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All  current  low-earth-orbit  satellites  are  simply  placed 
in  orbit  and  then  tracked  as  their  orbits  decay  deeper  into 
the  atmosphere.  This  orbital  decay  is  due  to  forces  on  the 
satellite  produced  by  atmospheric  drag  and  the  non- 
spherical  earth. 

Modified  exponential -type  models  for  the  atmosphere 
produce  sufficiently  accurate  drag  predictions  to  track  the 
satellites. 

New  concepts  in  scientific  and  military  uses  of  space 
require  greater  accuracy  in  predicting  atmospheric  drag 
effects  then  produced  by  the  modified  exponential  models. 
These  new  uses  include  the  Aero-assisted  orbital  Transfer 
Vehicle  (AOTV) ,  orbiting  space  stations,  the  Trans 
Atmospheric  Vehicle  (TAV) ,  and  the  strategic  defense 
initiative  satellite  constellation. 

The  environment  of  all  of  these  vehicles  is  called  the 
thermosphere,  an  atmospheric  layer  studied  in  depth  only 
recently.  An  accurate  empirical  model  of  the  thermosphere 
has  also  been  produced. 

One  idea  under  investigation  at  the  Naval  Surface 
Weapons  Center  is  the  placement  of  many  low-earth-orbital 
defensive  satellites  in  precisely  controlled  and  maintained 
orbits.  This  set  of  defensive  satellites  is  termed  a 
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constellation.  All  low  earth-orbital-satellites. drift  from 
their  initial  orbits  due  to  perturbating  accelerations 
caused  primarily  by  the  aspherical  earth  and  atmospheric 
drag.  In  order  to  be  functional,  the  orbits  of  the 
satellites  in  an  SDI  constellation  must  be  exactly 
maintained.  The  drift  of  a  satellite  must  be  cancelled  out 
by  using  an  onboard  propulsion  system  to  produce  forces  to 
offset  the  perturbating  accelerations. 

One  model,  proposed  by  Dr.  Dan  Parks  of  the  Naval 
Surface  Weapons  Center,  calculates  the  mission  lifetime 
expected  of  a  satellite  doing  intrack  micro-thrusting  to 
overcome  atmospheric  drag  effects.  The  end  of  mission  life 
corresponds  to  depletion  of  all  the  maneuvering  fuel.  This 
model  is  here-after  referred  to  as  the  propellant  longevity 
model.  The  propellant  longevity  model  incorporates  a  Bessel 
function  expansion  for  the  fuel-mass  decrement  on  each 
revolution.  However,  the  model  does  not  calculate  the  fuel 
necessary  to  cancel  out  orbital  drift  due  to  the  aspherical 
earth . 

This  thesis  gives  a  verification  of  the  low-altitude 
earth  satellite  propellant  longevity  prediction  model 
developed  by  Dr.  Parks.  The  theory  for  the  model 
development  itself  is  presented  in  Naval  Surface  Weapons 
Center  Technical  Report  (NSWC  TR)  83-243  (Parks  [Ref.  13). 
Computer  code  for  the  model  was  developed  at  NSWC  under  the 
direction  of  Dr.  Parks.  In  addition  to  the  computer  code  of 
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the  model,  several  other  NSWC  publications  with  applications 
to  the  model  were  forwarded  to  the  Naval  Postgraduate 
School . 

The  Thesis  is  organized  as  follows.  Chapter  I  presents 
a  general  introduction  to  the  propellant  longevity  model  and 
discusses  the  changes  made  to  the  computer  code  for  use  at 
the  Naval  Postgraduate  School.  A  basic  description  of  the 
orbital  mechanics  necessary  to  operate  the  model  is  provided 
in  Chapter  II.  Chapter  II  also  presents  an  overview  of  the 
atmospheric  physics  necessary  to  understand  the  density  and 
compositional  variations  in  the  thermospheric  layer  of  the 
thermosphere.  A  description  of  the  propellant  longevity 
model  is  given  in  Chapter  III.  Instructions  for  operation  of 
the  computer  code  are  also  included  there.  Chapter  IV 
contains  the  results  of  the  verification  process  sensitivity 
analysis  and  reasonableness  tests.  A  summary  is  given  in 
Chapter  V.  Appendix  A  contains  the  propellant  longevity 
Fortran  computer  code  for  the  IBM  3033.  Appendix  B  shows 
1  input  sets  and  their  corresponding  outputs.  Graphical 

presentation  of  the  tabular  results  is  included. 

The  coding  of  the  model  by  NSWC  was  completed  in  Fortran 
for  the  Control  Data  Corporation  (CDC)  6700  computer. 
Because  of  software  and  hardware  differences,  extensive 
modifications  to  the  model  coding  were  required  to  install 
the  program  as  an  operational  model  on  the  IBM  3033  computer 
at  the  Naval  Postgraduate  School.  To  optimize  the  utility 
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of  the  model,  some  adjustments  were  also  made  to  the 
original  computer  coding  of  the  propellant  longevity  model. 
Great  care  was  taken  when  modifying  the  code  to  preserve  the 
model  theory  as  contained  in  NSWC  TR  83-243. 

The  alterations  to  the  computer  program  fall  into  two 
categories:  changes  to  transport  and  implement  the  model  on 
our  IBM  3033,  and  changes  made  to  mainstream  the  model 
verification  process.  Changes  to  transport  the  model  from 
the  CDC  environment  to  the  IBM  environment  involved  some 
logical  revision  and  syntax,  hardware  word-size 
considerations . 

The  hardware  word-size  problems  were  extensive.  The  CDC 
6700  uses  a  60-bit  word  size  for  single  precision 
calculations,  whereas  IBM’s  standard  is  a  30-bit  word.  The 
propellant  longevity  model  requires  60  bit  precision.  Thus 
all  real  numbers,  library  functions,  and  input  constants 
were  promoted  to  double  precision  accuracy.  Another 
precision-related  problem  involved  the  use  of  fractional 
exponentiation.  All  one-half  exponents  were  changed  to  call 
the  SQRT  function  to  increase  precision  and  speed. 

Syntactical  changes  were  mostly  related  to  the  software 
environment  transportation  from  CDC  Fortran  to  IBM  Fortran 
77.  Library  functions  calls,  and  input  and  output  format 
statement  changes,  were  the  bulk  of  the  syntactical  changes 
made.  (Examples  of  library  function  changes  are  the  ARQTAN 
function  changed  to  the  AT AN 2  function,  the  BESI  Bessel 


function  call  changed  to  the  MMBSIR  Bessel  function  call  and 
so  forth) . 

Examples  of  changes  in  the  logic  to  in  order  to 
transport  and  implement  the  model  code  include  using  the 
Fortran  77  block  IF-THEN -ELSE  statement  for  the  CDC  Fortran 
IF  statement  and  the  use  of  DIMENSION  statements  in  each 
subroutine  that  is  processing  a  particular  array. 

The  second  category  of  changes  to  the  program  addressed 
mainstreaming  the  model  verification  process.  The  model 
code  as  down-loaded  from  the  NSWC  tape  contained  little 
code- internal  documentation.  Looking  towards  future 
implementation  and  maintenance  of  the  model  code,  extensive 
documentation  was  added  as  code  internal  comments.  Most  of 
this  documentation  was  based  on  the  research  literature  of 
NSWC.  Extensive  reference  to  the  current  publications  is 
made  in  the  code  documentation  itself. 

Two  changes  were  made  in  the  way  the  code  solves  the 
model  in  order  to  implement  the  theory  with  greater 
precision.  The  first  change  turned  on  the  subroutine  THRST 
which  was  designed  to  allow  for  preplanned  changes  to  the 
semi-major  axis  of  the  orbit  at  the  expense  of  onboard 
maneuvering  fuel.  (The  subroutine  code  down-loaded  from 
tape  did  not  return  the  fuel  consumed  for  the  orbit  adjust, 
record  the  resultant  change  in  the  semi-major  axis,  or 
receive  the  correct  intended  adjust  parameter.  All  these 
functions  are  performed  in  the  IBM  3033  version.) 


A  second  change  in  the  code  alters  the  way  in  which  the 
satellite  velocity  relative  to  the  atmosphere  is  computed. 
In  the  original  code  this  parameter  had  to  be  entered  as 
part  of  the  input  file  and  was  assumed  to  be  constant. 
Since  this  relative  velocity  is  actually  a  function  of  other 
transient  orbital  parameters,  recomputation  is  now 
automatically  accomplished  for  each  satellite  revolution 
based  on  current  values  of  the  changing  orbital  parameters. 
In  place  of  the  satellite-to-atmospheric-relative-velocity 
input,  the  user  now  inputs  the  atmospheric-to-earth  relative 
rotation  rates.  The  physics  associated  with  this  change  is 
discussed  later  on,  in  the  section  describing  the 
atmospheric  model. 

A  few  other  changes  to  the  propellant  longevity  code 
customize  the  output,  or  state  the  input  element  initial 
conditions.  Several  control  inputs  now  allow  for  a  range  of 
model  predictions.  These  are  explained  in  Chapter  III. 

The  actual  model  verification  is  accomplished  in  three 
steps.  First  does  the  model  address  the  original  problem? 
Second,  does  the  model  meet  the  criteria  of  common  sense? 
Finally,  how  well  does  the  model  solution  compare  when 
tested  against  real-world  data? 

No  real-world  data  sets  exist  for  testing  the  propellant 
longevity  model.  The  only  low-earth-orbit  satellite  having 
a  propulsion  system  to  maintain  orbital  parameters  is  highly 
classified  under  the  auspices  of  the  United  States  Air 


Force.  The  data  from  this  program  is  not  available  in  an 
unclassified  state.  In  lieu  of  actual  data,  several 
additions  to  the  program  code  were  made  to  assist  in  the 
model  verification  process.  The  major  adjustment  is  the 
addition  of  a  feature  to  allow  for  a  sensitivity  analysis. 
The  subroutine  SNALYS,  allows  the  user  to  linearly  scan  any 
combination  of  the  20  input  elements,  from  a  user  specified 
minimum  to  any  maximum,  with  user  specified  granularity. 
The  output  of  the  sensitivity  analysis  can  be  descriptive 
and  tabular,  or  graphical  in  nature.  The  SNALYS  subroutine 
allows  for  a  plot  of  model  prediction  against  the  entire 
range  of  the  input  elements.  Results  of  input  element 
sensitivity  analysis  are  presented  in  the  Chapter  IV.  The 
sensitivity  analysis  feature  allows  for  optimum  orbital- 
element  and  satellite-ballistic-coefficient  design. 

Thought  was  given  to  possible  model  installation  on 
other  computers.  The  model  is  large  and  requires  about  20 
minutes  of  dedicated  central  processor  time  to  complete  an 
average  run.  As  such  the  model  is  rather  limited  to 
mainframe  computers.  Initial  calculations  based  on 
processor  time  indicate  an  IBM  AT  microcomputer  with  the 
requisite  80287  co-processor  and  additional  RAM  memory  might 
require  about  ten  hours  to  complete  an  average  run. 

This  chapter  presented  an  overview  and  motivation  for 
the  propellant  longevity  model.  Changes  made  to  the 
computer  program  at  the  Naval  Postgraduate  School  were  also 


This  chapter  presents  the  orbital  mechanics  and 
atmospheric  concepts  required  for  the  understanding  and 
operation  of  the  propellant  longevity  model.  The  orbital 
mechanics  presented  in  the  first  section  discusses  only  the 
orbital  element  sets,  concepts,  and  coordinate  systems 
specific  to  the  model.  In  the  second  section  of  this 
chapter.  Thermosphere  specific  atmospheric  concepts  and 
cycles  are  presented. 

A.  ORBITAL  MECHANICS 

This  section  covers  the  orbital  mechanics  for  operating 
the  mass  decrement  model.  If  the  coordinate  and  orbital 
element  systems  discussed  are  familiar,  this  basic  orbital 
mechanics  section  can  be  skipped.  Orbital  mechanics 
references  are  cited  in  the  text. 

Early  celestial  mechanics  models  were  based  on 
increasingly  accurate  observations.  As  observations  and 
record  keeping  improved  so  did  the  empirical  model  fit  to 
the  observations.  The  model  that  fathered  modern  orbital 
mechanics  was  the  empirical  model  derived  by  Johannes 
Kepler,  1571-1630.  In  1619  Kepler  formulated  three  laws 
(models)  of  planetary  motion.  These  models  were  based  on 
life  long  observations  recorded  by  Tycho  Brahe,  1564-1601. 
Kepler's  models  were  formulated  as  follows: 


1.  The  orbit  of  each  planet  is  an  ellipse  with  the  sun  at 
one  focus. 


2.  The  line  joining  a  planet  to  the  sun  sweeps  out  equal 
areas  in  equal  times. 


3.  The  square  of  the  period  of  a  planet  is  proportional 
to  the  cube  of  its  mean  distance  from  the  sun. 


Later,  Sir  Isaac  Newton,  1643-1727  validated  Kepler's 


laws  by  establishing  their  equivalence  with  his  law  of 


universal  gravitational  attraction. 


Newton's  model  for  gravitational  attraction  states  that 


two  bodies  attract  with  a  force  that  is  proportional  to  the 


product  of  their  masses  and  inversely  proportional  to  the 


square  of  the  distance  between  them.  Symbolically,  the 


force  on  the  mass  m  due  to  the  mass  M  is  given  by: 


Eg  =  - 


GMm  r_ 
2  r 


where  G  is  the  universal  gravitational  constant  and  r  is  the 


vector  from  the  center  of  mass  of  M  to  the  center  of  mass  m. 


r  is  the  magnitude  of  the  distance  vector  r. 


A  body  experiences  an  acceleration  due  to  the  resultant 


of  all  the  forces  acting  on  the  body.  If  only  a 


gravitational  force  acts  on  the  body,  then  the  acceleration 


of  the  body  is  given  by: 


G(M+m)  ~ 

- 3  r 

r 
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In  this  model,  r"  is  the  acceleration  vector  on  m  caused  by 
gravitational  attraction  of  M  on  m.  The  remaining  variables 
are  defined  as  in  the  previous  model.  If  M  is  much  greater 
then  m  (as  in  an  earth  satellite  system),  G(M+m)  is  very 
nearly  GM.  Another  constant  can  then  be  defined  called  the 
gravitational  potential  parameter  u: 

u  =  Gm  ( 3 ) 

The  gravitational  acceleration  of  a  satellite  in  earth 
orbit  can  then  be  expressed  as: 

r"  +  u/r3  *  r  -  0  (4) 

The  gravitational  field  is  conservative.  This  means 
that  an  energy  constant  of  motion  can  be  defined  for  a  body 
moving  in  a  gravitational  field.  For  a  two  body,  an- 
atmospheric,  elliptic  orbit  system  the  specific  mechanical 
energy  is: 

£  -  vt2/2  -  u/r  (5) 

In  this  model  e,  is  the  specific  mechanical  energy,  vt  is 
the  satellite  tangential  speed,  and  u  and  r  are  as 
previously  defined.  Note  that  E  is  the  sum  of  the  kinetic 
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energy  per  unit  mass  of  the  satellite  and  its  potential 
energy  per  unit  mass. 

The  model,  (5),  requires  satellites  at  specific  radii  to 
have  specific  tangential  speeds.  As  r  increases,  v  must 
decrease  to  maintain  £  constant.  This  orbital  energy  model 
is  important  when  modelling  the  effects  of  atmospheric  drag 
which  is  the  major  topic  of  the  next  section. 

Newton's  inverse  square  law  of  gravitational  force 
requires  gravitational  field  orbits  to  be  elliptical  in  the 
shape.  The  mathematics  of  the  ellipse  is  important  in  the 
description  of  orbits  and  Figure  1  presents  some  of  the 
important  parameters. 

The  parameters  shown  in  Figure  1  are  defined  as  follows: 
a:  the  semi-major  axis, 
b:  the  semi-minor  axis, 
e:  the  eccentricity  where 

e  =  SQRT  [1  -  (b/a)2]. 

C:  the  ellipse  center. 

O:  the  ellipse  focus  (earth  center) . 
p:  the  semi-latus  rectum. 

r:  the  distance  from  a  position  on  the  ellipse 

the  focus  0. 

v:  the  true  anomaly. 

rp:  the  radius  at  perigee  (in  an  earth  system) . 
ra:  the  radius  at  apogee  (in  an  earth  system). 
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Figure  1.  The  Elliptical  Orbit 

A:  the  apogee  (point  farthest  from  the  focus) . 

£:  the  perigee  (point  closest  to  the  focus. 

No  real  orbits  occur  in  just  two  dimensions.  To  describe 
real  orbits  a  third  dimension  must  be  added.  A  suitable 
reference  must  also  be  chosen.  One  of  the  most  convenient 
reference  frames  is  the  geocentric  reference  depicted  in 
Figure  2.  In  this  system  a  fixed  nonrotating  earth  forms 
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-  the  center  of  the  celestial  sphere.  The  earth's  projected 
north  pole,  south  pole,  and  equator  become  the  north 
celestial  pole  (NCP) ,  the  south  celestial  pole  (SCP) ,  and 
celestial  equator  respectively. 


Figure  2.  The  Geocentric  Celestial  Sphere 

To  account  for  the  apparent  motion  of  the  celestial 
sphere  as  viewed  from  the  earth,  a  reference  point  on  the 
celestial  sphere  must  be  chosen.  Traditionally  this  has 
been  the  First  Point  of  Aries,  describing  the  constellation 
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where  the  vernal  equinox  occurred  when  the  celestial 


reference  was  defined.  The  vernal  equinox  does  not  now 
occur  in  the  constellation  Aries.  Tradition  still  uses  the 
first  point  of  Aries  as  the  Celestial  Meridian  (CM)  of 


reference. 


Figure  3  illustrates  the  vernal  equinox 


reference.  The  vernal  equinox  is  defined  as  the  time  when 
the  sun  is  directly  over  the  equator  at  noon  in  the  spring. 
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Figure  3.  The  Vernal  Equinox 
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Figure  4  is  a  representation  of  an  orbit  in  real  space 
using  the  geocentric  system.  The  points  where  the 
projection  of  the  satellites  orbit  on  the  celestial  sphere 
crosses  the  celestial  equator  are  called  the  nodes  of  the 
orbit.  The  ascending  node  (N)  is  where  the  projected  orbit 
crosses  the  celestial  equator  heading  for  the  north 
celestial  pole.  The  descending  node  is  where  the  projected 
orbit  crosses  the  celestial  equator  heading  for  the  south 
celestial  pole.  The  angular  measurement  from  the  reference 
celestial  meridian  to  the  ascending  node  along  the  celestial 
equator  is  called  the  longitude  of  the  ascending  node,  and 
usually  given  the  symbol  .  The  angular  measurement  from 
the  ascending  node  to  the  point  of  perigee  (P)  along  the 
projected  orbit  is  called  the  argument  of  perigee, 
symbolized  by  w .  Both  the  longitude  of  the  ascending  node 
and  the  argument  of  perigee  range  from  0  to  360  degrees. 
The  angle  from  the  plane  of  the  celestial  equator  to  the 
plane  of  the  projected  orbit  is  called  the  inclination  and 
denoted  by  i.  Values  for  the  inclination  range  between  0 
and  180  degrees. 

Inclinations  of  0  and  180  degrees  are  termed  equatorial. 
In  equatorial  orbits  the  longitude  of  the  ascending  node  is 
not  defined.  Orbits  with  an  inclination  of  90  degrees  are 


called  polar.  Prograde  orbits  are  those  with  an  inclination 
of  between  0  and  90  degrees  and  retrograde  orbits  have  an 
inclination  range  of  90  to  180  degrees. 
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Figure  4.  A  Geocentric  Orbital  System 


Historically,  it  was  difficult  to  tell  much  about  the 
orbit  specific  parameters  from  the  earth  vantage  point. 
However,  a  value  that  could  be  measured  was  the  time  for  one 
complete  orbit.  If  the  radius  vector  moved  uniformly  along 
the  projected  orbit  with  the  average  rate  of  mean  motion  n 
(n  is  also  referred  to  as  the  mean  angular  velocity)  at  an 

epoch  T,  it  would  be  at  a  position  M  called  the  mean 
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anomaly.  The  mean  anomaly  has  a  range  of  0  to  360  degrees. 
It  is  measured  from  the  perigee  along  the  projected  orbit. 
The  mean  anomaly  corresponds  to  positions  on  the  orbital 
ellipse.  An  eccentric  anomaly  (E) ,  is  defined  as  the  angle 
between  the  projected  position  and  the  point  of  perigee 
subtending  at  the  center  of  the  ellipse.  Figure  5 
illustrates  eccentric  anomaly. 


Figure  5.  Orbital  Anomalies 
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The  difference  between  the  eccentric  anomaly  and  the  mean 
anomaly  is  important.  The  mean  anomaly  is  a  mathematical 
concept  corresponding  to  an  imaginary  position  of  a 
satellite  past  the  perigee  if  it  moved  at  an  average  angular 
velocity.  The  eccentric  anomaly  is  the  actual  projected 
angular  position  of  a  real  satellite  on  the  celestial 
sphere.  The  eccentric  and  mean  anomalies  are  equal  in 
circular  orbits  (e  =  0)  .  The  eccentric  and  mean  anomalies 
are  related  by  Kepler's  equation: 


E  -  e*sin  E 


The  solution  to  model  (6)  has  been  the  source  of  much 
investigation.  By  a  solution  it  is  meant  for  a  given  M, 
determine  e  and  E.  In  1817  Friedrich  Wilhelm  Bessel,  1784- 
1846  invented  Bessel  functions  to  provide  an  infinite  series 
approximation  for  the  solution  of  the  model  (6)  .  Newton 
developed  a  converging  algorithm  for  its  solution.  The  mass 
decrement  submodel  of  the  propellant  longevity  model  uses 
Newton's  method  for  the  solution  of  Kepler's  equation. 

All  anomalies  are  concerned  with  the  measurement  of  the 
position  of  a  body  in  orbit.  All  anomalies  are  measured 
angularly  from  the  point  of  perigee  in  the  orbital  plane, 
and  vanish  at  perigee. 

When  e  and  E  have  been  found  from  Kepler's  equation  (6), 
then  using  the  equations  for  an  ellipse  in  Taff  [Ref.  2:p. 
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67],  the  remaining  elliptical  parameters  can  be  solved  for 
when  any  one  of  a,  b,  rp,  or  ra  is  known. 

The  minimum  number  of  parameters  needed  to  describe  an 
orbit  is  called  an  orbital  element  set.  Five  parameters  are 
needed  to  describe  an  orbit.  One  more  is  needed  to  position 
an  object  in  the  orbit.  Many  orbital  element  sets  have  been 
used  and  are  in  use.  The  classical  orbital  elements  are: 
semi-major  axis  (a) 
eccentricity  (e) 
inclination  (i) 

longitude  of  the  ascending  node  (ft) 

argument  of  perigee  (u) 

time  of  perigee  passage  or  epoch  (T) . 

The  Keplerian  set  replaces  the  time  of  perigee  passage  with 
the  mean  anomaly  (M)  .  The  orbital  element  set  used  in  the 
mass  decrement  submodel  is  the  Kozai  mean  element  set,  very 
similar  to  the  Keplerian  set.  The  name  "Kozai  mean  element 
set"  refers  to  the  Kozai  perturbation  approximations  used  to 
solve  for  the  Keplerian  orbital  element  set.  The 
perturbation  theory  approximation  of  Y.  Kozai  is  discussed 
later. 

Up  to  now  all  orbital  considerations  have  been  limited 
to  the  classical  two-body  problem.  The  two  body  problem 
limits  the  forces  on  the  orbiting  body  to  the  gravitational 
force.  The  two-body  problem  also  requires  all  masses  be 
concentrated  at  a  single  point.  A  real  earth  satellite 
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experiences  many  forces  resulting  in  many  discrete 
accelerations.  The  theory  of  perturbations  was  developed  to 
account  for  the  motion  of  real  satellites. 

In  perturbation  theory  the  original  orbit,  with  only  the 
point -mass  gravitational  acceleration,  is  called  the 
osculating  orbit.  Orbital  elements  describing  the 
osculating  orbit  are  the  osculating  element  set.  The 
osculating  orbit  is  then  adjusted  periodically  to  account 
for  the  action  of  all  the  nonpoint-mass  gravitational 
accelerations . 

Accelerations  result  from  forces.  Examples  of 
perturbating  forces  on  an  earth  satellite  are:  oblate  earth 
gravitational,  atmospheric  drag,  sun  gravitational,  moon 
gravitational,  solar  radiation  pressure,  and  satellite  motor 
thrust . 

To  remain  in  a  closed  orbit,  the  point-mass 
gravitational  force  must  be  the  dominant  force  on  the 
satellite.  For  low  earth  orbit  satellites,  the  oblate  earth 
gravitational  force  is  the  next  greatest  force  and  is  termed 
the  dominant  perturbation.  Satellites  below  600  kilometers 
in  altitude  experience  an  atmospheric  drag  force  of  the  same 
order  of  magnitude  as  the  oblate  earth  force.  The 
atmospheric  drag  and  oblate  earth  effects  are  so  large  for 
low  earth  orbit  satellites  that  the  effects  of  the  other 
perturbations  (except  motor  thrust)  are  lost  in  the 
background  error. 


25 


Central  to  the  problem  of  perturbation  theory  is  the 
that  the  three  body  problem  has  never  been  solved 
mathematically  in  general  closed  form.  Most  modern 
approximations  involve  power  series  or  Legendre  polynomial 
expansions.  The  state  of  the  art  in  orbit  perturbation 
prediction  only  uses  only  the  first-order  terms  of  these 
expansions,  termed  first-order  effects.  An  exception  is  a 
theory  advanced  in  1959  by  Y.  Kozai.  His  approximations  are 
computationally  efficient  and  allow  for  the  use  of  more  of 
the  oblate  earth  effects  (called  J2,  J3,  and  J4  terms  and 
discussed  later).  Kozai' s  theory  is  flawed  in  the  higher- 
order  terms  due  to  the  introduction  of  mathematical 
singularities.  Nevertheless,  proponents  for  its  use  argue 
that  it  does  work  and  predicts  the  oblate  earth  effects  to  a 
reasonable  degree  of  accuracy.  The  mass  decrement  submodel 
of  the  propellant  longevity  model  presented  in  this  thesis 
uses  Kozai's  theory. 

The  gravitational  potential  parameter  u  was  presented  as 
singularly  determining  the  gravitational  potential  in  the 
earlier  development  of  the  gravitational  acceleration  model. 
In  earth  systems,  the  gravitational  potential  is  a  complex 
function  of  the  earth's  mass  distribution. 

The  earth  is  not  a  point  mass.  Neither  is  it  a 
homogenous  sphere.  It's  true  mass  distribution  is  complex. 
The  real  distribution  of  matter  on  the  earth  forms  a  body 


termed  the  geoid.  This  geoid  determines  exactly  the 
gravitational  potential  function. 

The  determination  of  the  geoid  is  based  on  ground  and 
satellite  measurements.  Many  military  programs  depend  on 
the  precise  position  of  a  point  on  the  surface  of  the  earth. 
The  military  has  led  the  way  in  geoid  determination  in  the 
past  two  decades. 

Geoidal  studies  reveal  important  facts  about  the 
aspherical  shape  of  the  earth.  A  fluid  spinning  sphere 
should  exhibit  a  bulge  at  the  equator  and  flattening  at  the 
poles.  This  is  termed  oblateness.  The  oblate  deviation  is 
the  dominant  deviation  for  the  earth  as  a  sphere. 
Measurements  of  the  earth's  oblateness  indicate  that  it  is 
more  oblate  than  the  current  spin  rate  could  produce  by 
itself.  This  excess  oblateness  indicates  two  factors:  the 
earth's  core  is  more  solid  than  previously  thought,  and  the 
earth's  spin  has  slowed. 

The  earth  deviates  in  other  ways  from  being  a  sphere. 
Since  all  deviations  can  be  expressed  in  terms  of  sines  and 
cosines  they  are  referred  to  as  spherical  harmonics. 
Harmonics  are  divided  into  categories.  Zonal  harmonics  are 
symmetric  about  the  earth's  spin  axis  and  are  not  longitude 
dependent.  Sectorial  harmonics  are  only  longitude 
dependent.  Harmonics  dependent  on  both  longitude  and 
latitude  are  termed  tesseral.  Odd  numbered  harmonics  are 
antisymmetric  about  the  equatorial  plane;  even  numbered 
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harmonics  are  symmetric  about  the  equatorial  plane. 
Harmonics  are  also  classified  by  number.  The  number  refers 
to  the  subscript  of  the  J  coefficient  (defined  later) . 

Bate  et  al.,  [Ref.  3],  offers  a  simplified  development 
of  the  gravitational  potential  function.  If  the 
gravitational  potential  function  <J>  of  a  geoid  is  known,  the 
acceleration  due  to  this  gravitational  function  is  modeled 
by: 

r"  =  (7) 

One  common  expression  for  the  gravitational  potential  or 
geoid  function  <j>  is: 

*  -  7II  -  I  J  (~) n  Pn  sin  L]  (8) 
r  n=2  n  r  n 

In  this  model  u  is  the  gravitational  potential  parameter,  Jn 
is  an  experimentally  determined  J  coefficient,  re  is  the 
radius  of  the  earth,  Pn  is  a  Legendre  polynomial  of  order  n, 
L  is  the  geocentric  latitude,  and  r  is  the  radius  of  the 
satellites  position.  Taking  the  gradient  of  the  potential 
function  yields  an  expression  for  the  acceleration  due  to 
the  gravitational  potential.  Bate  et  al.,  [Ref.  3:pp.  419- 
423]  details  the  resultant  potential  induced  acceleration  to 
the  first  seven  terms. 
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Recent  studies  have  determined  the  geoid  function  out 
through  the  J44  term.  Values  used  by  the  mass  decrement 
model  by  the  subroutine  GEOP  are  1082 . 76*10“6,  -2.55*10”6, 
and  1 . 56*10”® ,  for  the  32,  32,  and  J4  terms  respectively. 
The  32  and  J4  terms  are  even  numbered,  and  refer  to  even 
harmonics  symmetric  about  the  equator.  The  32  term  is  odd 
numbered,  and  refers  to  odd  harmonics  which  are 
antisymmetric  about  the  equator.  The  32,  32,  and  J4  terms 
are  all  zonal. 

Current  perturbation  theory  expresses  $  to  the  32 
term.  The  theory  advanced  by  Kozai  in  1959  appears  to  use 
the  32,  32,  and  J4  coefficients  in  a  first-order 
perturbation  theory  approximation  (because  singularities 
occur  in  terms  higher-order  than  32) .  Taff  [Ref.  2:pp.  332- 
342]  takes  strong  exception  to  the  theory. 

The  subroutine  GEOP  in  the  propellant  longevity  model 
uses  Kozai' s  theory  to  calculate  the  geoidal  effects  on 
orbital  parameters  through  the  first  three  zonal 
coefficients  (J2,  J3,  and  J4  terms).  Their  effect  is 
calculated  once  each  orbital  revolution,  and  the  Kozai  mean 
element  set  is  adjusted  accordingly.  While  using  the  model, 
it  is  important  to  understand  that  as  the  altitude  of  the 
satellite  decreases,  The  aspherical  earth  perturbations 
become  much  greater.  The  orbital  element  drift  due  to  the 
aspherical  earth  is  thus  much  faster  in  lower  altitude 
orbits. 
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In  summary,  the  computer  program  of  the  propellant 
longevity  model  uses  three  orbital  element  sets  in  a 
geocentric  reference.  The  input  orbital  element  set  is  the 
Brouwer  mean  element  set  of  the  Naval  Space  Surveillance 
System  (NAVSPASUR) 1  system.  Chapter  III  and  Appendix  A 
explain  the  input  elements  in  detail  in  the  input 
description.  The  Brouwer  mean  elements  are  first  converted 
into  a  Brouwer  osculating  set.  These  osculating  orbital 
elements  are  then  transformed  into  a  Keplerian-like  Kozai 
mean  set.  Half  of  the  computer  program  of  the  propellant 
longevity  model  transforms  the  Brouwer  mean  element  set  into 
the  Kozai  mean  set.  These  element  set  transformations  are 
used  to  prepare  real-world  orbital  element  sets  for  use  in 
the  propellant  longevity  predictions.  Once  the  computer 
program  has  solved  for  the  Kozai  mean  element  set  all 
perturbation  effects  result  in  adjustments  to  the  Kozai  mean 
set.  Adjustments  are  made  once  per  revolution. 

B.  THE  ATMOSPHERE  MODEL 

Central  to  the  validation  of  any  model  is  an 
understanding  of  the  system  to  be  modeled.  The  model  under 
investigation  determines  the  propellant  longevity  for  low 
earth  orbital  satellites  doing  in  track  micro-thrusting  to 

^■The  NAVSPASUR  system  publishes  a  semiannual  report  on 
the  1000+  orbiting  objects  they  track  in  The  Satellite 
Situation  Summary.  The  orbits  are  described  in  a  five  card, 
80  column,  Fortran  style  input  format.  The  report  is 
available  on  request  from  the  Commanding  Officer,  Naval 
Space  Surveillance  System,  Dahlgren,  Virginia. 
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overcome  the  effects  of  atmospheric  drag.  The  fuel  required 
by  a  particular  rocket  motor  on  a  satellite  to  perform  a 
known  amount  of  work  can  be  computed  from  the  rocket 
equation.  The  problem  occurs  in  computing  the  amount  of 
work  required. 

To  compute  the  required  work  by  the  satellite  propulsion 
system  to  overcome  atmospheric  drag  requires  knowledge  of 
the  atmospheric  environment  of  the  satellite.  For  low  earth 
satellites  the  atmospheric  environment  is  the  thermosphere. 
The  physics  of  the  thermosphere  is  discussed  in  this 
section.  Before  discussing  the  thermosphere,  some 
preliminary  thermospheric  properties  are  presented. 

For  low  earth  orbital  satellites,  the  drag-induced 
change  in  the  orbital  elements  is  significant.  The  oblate 
earth  (referred  to  as  the  "J2  term")  orbital  perturbation  is 
normally  the  dominant  perturbation  term.  The  perturbations 
due  to  drag  are  of  the  same  order  of  magnitude  as  those  due 
to  the  oblate  earth  for  orbits  below  600  kilometers. 

The  drag  force  caused  by  the  satellite  moving  through 
the  atmosphere  always  acts  against  the  motion  of  the 
satellite.  Neglecting  the  lift  effects  due  to  the 
satellite's  shape  and  attitude,  the  drag  force  causes  work 


to  be  done  on  the  atmosphere  by  the  satellite. 

*5rora  Bate  e  al.,  [Ref.  3:pp.  15-16]  a  satellite  in  orbit 
has  a  constant  specific  mechanical  energy  £  given  by 
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where  £  is  the  mechanical  energy  per  unit  mass,  vt  the 
satellite  speed,  u  the  gravitational  potential  constant,  and 
r  the  radius  of  the  orbit. 

The  above  model  says  that  the  mechanical  energy  per  unit 
mass  of  a  body  in  an  an-atmospheric,  two  body,  circular 
orbit  is  the  sum  of  its  kinetic  energy  per  unit  mass  and 
gravitational  potential  energy  per  unit  mass. 

Conservation  of  energy  requires  that  the  work  done  on 
the  atmosphere  by  the  satellite  must  reduce  the  orbital 
energy  of  the  satellite.  The  differential  element  of  work 
done  on  the  atmosphere  by  the  satellite  results  in  a 
differential  reduction  in  the  specific  kinetic  energy  of  the 
satellite.  Assuming  the  specific  mechanical  energy  is  a 
constant,  the  radius,  r,  of  the  orbit  must  therefore 
decrease.  As  the  radius  decreases,  the  velocity  tends  to 
increase  in  order  to  maintain  the  an-atmospheric 
conservation  of  specific  mechanical  energy. 

In  non-circular  (real  elliptical)  orbits,  the 
atmospheric  effects  are  greatest  at  perigee  (Figure  l)  . 
Orbital  mechanics  requires  a  differential  reduction  in 
orbital  velocity  at  perigee  to  result  in  a  reduction  in  the 
semi-major  axis  and  orbital  eccentricity.  The  major  effect 
of  transatmospheric  satellite  motion  is  to  reduce  the  semi¬ 
major  axis  and  the  eccentricity  of  the  elliptical  orbit. 
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Taff  [Ref.  2:pp.  352-353]  develops  the  mathematical 
equations  for  the  drag-induced  changes  in  the  semi-major 
axis,  eccentricity,  argument  of  perigee  (Figure  4) ,  and  mean 
anomaly  (Figure  5)  .  The  changes  in  the  argument  of  perigee 
and  the  mean  anomaly  are  shown  to  be  periodic.  The  net 
effect  on  the  orbit  due  to  atmospheric  drag  is  to 
circularize  the  orbit  closer  to  the  earth,  with  little 
change  to  the  other  orbital  elements. 

An  accepted  model  for  the  magnitude  of  the  drag  force  as 
presented  in  Taff  [Ref.  2:p.  352]  is: 


*D 


A  pCj. 
2m 


rel 


(10) 


In  this  equation,  FD  is  the  drag  force  magnitude,  A  the 
effective  satellite  cross-sectional  area,  CD  the  drag 
coefficient,  m  the  satellite  mass,  p  the  atmospheric 
density,  and  vrel  the  velocity  of  the  satellite  relative  to 
the  atmosphere. 

The  upper  atmosphere  is  dynamic.  The  changes  that  must 
be  modeled  to  achieve  a  valid  prediction  for  the  drag  force 
are  those  that  effect  the  drag  force  factors.  The  factors 
in  the  drag  force  model  that  are  dependent  on  the  atmosphere 
are  density,  a  drag  coefficient  (dependent  on  satellite 
physical  and  fluid-flow  characteristics) ,  and  the  relative 
velocity  of  the  satellite  in  the  atmosphere. 


To  model  the  atmospheric  effects  on  density,  drag 
coefficient,  and  relative  velocity,  some  background  in  upper 
atmosphere  physics  is  needed.  The  descriptive  presentation 
here  is  carefully  developed  to  be  relative  to  the  above 
three  drag  force  dependent  factors.  A  more  comprehensive 
treatment  of  atmospheric  physics  can  be  found  in  a  work  of 
Fleagle  [Ref.  4].  An  excellent  technical  treatment  of 
thermospheric  dynamics  can  be  found  in  the  works  of  Jacchia 
[Refs.  5,6,7]. 

Jacchia  formulated  the  most  widely  used  thermospheric 
model  currently  employed  by  professional  orbital  analysts 
(known  as  the  Jacchia  J60  model  or  simply  J60) .  His  latest 
comprehensive  empirical  thermospheric  model  (Jacchia  J77  or 
J77)  is  many  orbital  analysts'  preferred  refinement  for 
modeling  the  atmosphere. 

The  Jacchia  J77  model  is  highly  tabular,  uses 
considerable  computer  time,  and  requires  an  extensive  input 
file.  However,  these  negative  aspects  are  balanced  by  the 
orders-of-magnitude  increase  in  accuracy  of  density 
predictions.  Jacchia  [Ref.  7:p.  2]  presents  a  table  of 
residuals  for  predicted  versus  observed  densities  for 
various  satellites,  as  well  as  time  intervals  for  the 
Jacchia  J77  atmosphere  model.  All  density  residuals  are 
much  smaller  than  those  produced  by  exponential  models  in 
the  same  time  intervals.  When  extremely  accurate  density 
predictions  are  required  in  some  submodel,  the  use  of  the 


Jacchia  J77  atmosphere  model  may  well  be  worth  the 
additional  computer  time  and  expense. 

A  simple  model  of  the  atmosphere  considers  it  as  a 
homogeneous  mass  of  constant  density.  A  refinement  is  made 
by  modeling  the  density  as  decreasing  exponentially  with 
altitude.  The  exponential  model  predicts  pressure  and 
density  will  fall  one  tenth  for  every  ten  kilometer  increase 
in  altitude.  Prior  to  1960  this  exponential  model  was  used 
extensively  for  upper  air  calculations  involving  meteor 
trails.  With  the  advent  of  low  earth  orbital  satellites,  a 
greater  degree  of  accuracy  in  predicting  atmospheric  effects 
is  required.  The  predicted  satellite  decays  from  the 
exponential  model  are  very  different  from  the  observations. 
The  Jacchia  J60  atmosphere  model  was  formulated  to  improve 
decay  predictions.  It  was  based  on  measurements  form  high 
altitude  balloons,  sounding  rockets,  and  ground  stations. 

The  atmosphere  is  divided  roughly  into  concentric 
shells,  each  with  specific  properties.  The  properties  used 
in  this  subdivision  are  after  Fleagle.  [Ref.  4] 

From  the  early  days  of  atmospheric  studies,  temperature 
and  moisture  have  been  recorded.  The  division  system 
described  by  Fleagle  [Ref.  4]  is  based  primarily  on  a  height 
versus  temperature  profile. 

The  height  versus  temperature  profile  is  graphed  in 
Figure  6.  The  temperatures  and  altitudes  are  those  of  the 
1975  U.S.  standard  atmosphere.  This  temperature  behavior 
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mass.  This  region  is  characterized  by  near  linearly 
decreasing  temperature  with  altitude.  The  dynamics  are  due 
to  earth  surface  energy  transitions.  All  of  the  earth's 
weather  occurs  in  this  region.  The  height  of  the 
troposphere  is  variable  and  depends  on  both  latitude  and 
season.  Along  the  equator  the  maximum  altitude  is  16  to  18 
kilometers.  At  the  poles  the  troposphere  is  highly 
seasonal,  being  nearly  absent  during  the  winter  months  and 
extending  to  eight  to  ten  kilometers  in  the  summer.  Capping 
the  troposphere  is  a  region  of  constant  temperature  with 
altitude,  called  the  tropooause. 

Next  in  the  series  of  atmospheric  layers  is  the 
stratosphere .  The  stratosphere  is  characterized  by  linearly 
increasing  temperature  with  height.  The  physical  dynamics 
are  due  mainly  to  the  absorption  of  solar  ultra-violet 
radiation  by  the  ozone.  The  maximum  altitude  of  the 
stratosphere  is  about  50  kilometers  where  it  is  capped  by  an 
area  of  local  maximum  temperature  termed  the  stratopause. 

Above  the  stratopause  is  the  mesosphere.  This 
atmospheric  layer  exhibits  linearly  decreasing  temperature 
with  height.  The  energy  physics  of  the  mesosphere  are  not 
well  understood.  Bulk  mixing  of  the  gasses  still  dominates 
diffusion  as  thermodynamic  processes.  Any  molecular 


dissociation  of  atmospheric  molecules  produced  by  radiation 
absorption  is  rapidly  consumed  by  chemical  recombination  of 
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these  dissociated  species.  The  mesosphere  is  topped  by  the 
turbopause  at  an  altitude  of  about  90  kilometers. 

Above  the  turbopause  the  gross  physics  of  the  atmosphere 
changes  drastically.  This  region  is  the  environment  of  all 
current  low  earth  orbit  satellites.  Termed  the 
thermosphere .  it  extends  from  the  turbopause  to  an  altitude 
of  1000  to  2500  kilometers.  Here  bulk  molecular  mixing  is 
dominated  by  diffusion.  Lighter  gasses  diffuse  upwards  and 
heavier  ones  remain  at  the  lower  altitudes.  Disassociation 
of  molecular  species  becomes  increasingly  important,  this 
results  in  an  increase  of  individual  atomic  species. 
Standard  flow  equations  (Navier-Stokes)  are  extremely 
difficult  to  apply,  and  physical  descriptions  are  more 
discrete  (quantum  mechanical) .  Flight  by  flow-supplied  lift 
is  not  possible.  The  standard  methods  of  computing  drag 
coefficients  do  not  apply.  The  ideal  gas  approximation  is 
also  not  valid  (P*V  =  n*R*T  does  not  apply)  .  Above  the 
turbopause,  density  is  an  increasing  function  of 
temperature.  Only  about  one  millionth  of  the  atmosphere 
remains  above  the  tropopause. 

The  stratification  of  gaseous  components  in  the 
thermosphere  is  significant.  The  drag  coefficient  has  been 
shown  to  be  composition  dependent.  Thus,  as  the  gaseous 
composition  changes  the  drag  coefficient  must  be  adjusted  to 
reflect  the  observed  drag-induced  orbital  decays  (Jacchia 
[Ref.  7]).  Figure  7  shows  the  dominant  atmospheric  species 
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COMPOSITION  &  COEFFICIENT 

Figure  7.  Atmospheric  Composition 

by  altitude  with  the  accepted  drag  coefficients.  The  drag 
coefficients  are  predictions  obtained  analyzing  the  drag 
coefficient  needed  to  produce  the  observed  drag  effects  on 
orbital  decay  for  a  satellite  at  a  particular  altitude. 
Thus  the  drag  coefficient  is  a  proportionality  constant  in 
the  drag  force  model. 


For  altitudes  below  110  kilometers,  molecular  nitrogen 
is  the  main  atmospheric  gas.  Above  this  region  the 
dissociation  of  the  molecular  oxygen  by  ultra-violet  light 
absorption  causes  atomic  oxygen  to  become  increasingly  the 
main  gaseous  constituent.  At  an  altitude  of  500  kilometers, 
helium  replaces  atomic  oxygen  as  the  main  atmospheric  gas. 
At  altitudes  above  1500  kilometers  atomic  hydrogen  becomes 
the  dominant  gas  (see  Figure  7) . 

Traditional  research  has  assigned  a  drag  coefficient  of 
2.20  to  altitudes  between  200  and  400  kilometers.  As  the 
altitude  of  the  satellite  increases  above  500  kilometers, 
corresponding  to  satellite  movement  through  helium,  a  drag 
coefficient  of  at  least  3.00  must  be  assigned  to  maintain 
agreement  with  observation  (see  Figure  7) . 

Jacchia  [Ref.  7]  lists  the  observed  variational  changes 
in  density  and  composition  in  the  thermosphere.  They  are 
the  two  component  solar  activity  related,  semiannual, 
diurnal,  seasonal  latitudinal,  and  geomagnetic. 

The  two  component  solar  thermospheric  variations  are 
consequences  of  changes  in  solar  ultra-violet  (EUV)  flux. 
The  EUV  flux  changes  periodically  due  to  two  solar  cycles. 
These  are  the  27  day  rotational  or  disk  cycle  and  the  ll 
year  sunspot  or  active  area  cycle. 

EUV  flux  is  not  measurable  from  the  ground  due  to  almost 
complete  absorption  in  the  atmosphere.  Empirical  evidence 
has  suggested  a  strong  parallel  of  the  EUV  flux  and  the 


ground  measurable  10.7  centimeter  flux  (F10.7)  which  has 
been  used  in  some  successful  thermospheric  models.  The 
Jacchia  1960  model  used  the  10.7  centimeter  flux  as  an  input 
to  model  thermospheric  density  resulting  in  an  improvement 
of  decay  predictions  over  the  exponential  model.  (Indeed, 
it  is  only  recently  that  there  has  been  a  requirement  for 
greater  accuracy  in  drag  predictions.) 

The  thermosphere  responds  to  increased  EUV  flux  by 
absorbing  more  energy,  causing  increased  heating.  This 
increased  heating  causes  an  increase  in  temperature  and 
density.  The  thermospheric  response  to  changes  in  the  EUV 
flux  is  nonlinear.  Increases  in  EUV  flux  due  to  disk 
(rotational)  effects  produce  more  thermospheric  heating  then 
a  like  increase  in  flux  due  to  active  area  effects. 
However,  Jacchia  noticed  in  1973  that  the  disk  component  of 
the  EUV  flux  can  be  linearly  related  to  an  average  10.7 
centimeter  flux.  He  recommends  using  a  71  day  average 
corresponding  to  three  solar  rotations  [Ref.  7],  Tables  for 
the  10.7  centimeter  flux  have  been  developed.2  The  Jacchia 
1977  atmosphere  incorporates  the  average  10.7  centimeter 
flux  as  an  input,  used  as  an  indicator  of  the  thermospheric 
heating  capacity  of  the  sun's  radiation. 

20rr  [Ref.  8:Appendix  B]  lists  values  compiled  and 
predicted  by  the  NASA  Marshall  Space  Flight  Center.  Values 
range  from  a  minimum  of  7  3.3  to  a  maximum  of  2  4  2.9.  The 
units  for  the  10.7  centimeter  flux  are  10-22 
watts/meter2/Hz . 
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The  semiannual  temperature  and  density  variation  has 
been  modeled  in  many  ways.  Jacchia  [Ref.  7:pp.  47-50]  fits 
a  modified  harmonic  empirical  model  to  the  biannual, 
sinusoidal  cycle  of  observed  density.  Maxima  occur  in  April 
and  October.  The  Jacchia  1971  and  1977  models  for  the 
semiannual  density  variation  are  based  on  observations  of  12 
years  of  satellite  data.  The  observed  density  variations 
are  thought  to  be  a  consequence  of  the  earth-sun 
orientation.  The  mechanism  for  the  observed  effect  is 
unknown,  but  may  be  related  to  the  maximum  equatorial 
heating  at  the  equinoxes.  Equinoxes  (sun  at  maximum 
altitude  over  the  earth's  equator)  occur  in  late  March  and 
September.  Figure  8  illustrates  the  semiannual  bulge  in 
temperature  and  density  observed  just  after  the  equinoxes. 
Due  to  the  earth's  oblation,  the  bulk  mass  of  the  atmosphere 
is  concentrated  at  the  equator.  The  most  direct  exposure  of 
the  equatorial  atmosphere  occurs  at  the  equinoxes.  The 
temperature  and  density  thus  peak  shortly  after  periods  of 
maximum  heating. 

Seasonal-latitudinal  thermospheric  variations  are  again 
earth-sun  orientation  related.  Variations  in  density  and 
composition  are  dependent  on  the  season  and  latitude.  The 
Jacchia  1977  model  fits  modified  harmonic  empirical  models 
to  observed  data.  The  seasonal  latitudinal  variations 
interact  strongly  with  the  semiannual  variation,  complicat¬ 
ing  the  modeling.  The  winter  helium  bulge  is  the  best  known 
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Figure  8.  The  Semiannual  Atmospheric  Density  Bulge 

effect  of  this  variational  factor.  Figure  9  shows  the 
general  shape  change  of  the  temperature  and  density  bulge 
with  season.  Note  that  during  the  June  solstice  the  bulge 
is  centered  at  about  45  degrees  north  latitude  whereas 
during  the  December  solstice  the  bulge  is  centered  at  about 
45  degrees  south  latitude. 
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Figure  9.  The  Seasonal-Latitudinal  Atmospheric 
Density  Bulge 

The  diurnal  bulge  is  shown  in  Figure  10.  It  is  caused 
by  the  daytime  heating  of  the  side  of  the  atmosphere  facing 
the  sun.  Lagging  the  actual  (most  direct)  solar  exposure  by 
about  four  hours,  the  heated  side  shows  a  marked  increase  in 
temperature  and  density.  Many  current  models  track  and 
predict  the  position  of  the  diurnal  bulge. 
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Figure  10.  The  Diurnal  Atmospheric  Density  Bulge 


Geomagnetic  heating  of  the  thermosphere  is  not  well 


understood . 


At  times  of  high  geomagnetic  activity 


(corresponding  to  an  index  of  7)  large  heating  effects  are 
produced.  The  combined  effects  of  geomagnetic  and  EUV 
heating  can  produce  changes  in  density  of  three  orders  of 
magnitude  at  600  kilometers  (Jacchia  (Ref.  5:p.  82]). 
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In  summary,  the  thermosphere  is  a  very  dynamic  region  of 
the  atmosphere.  Many  of  the  observed  changes  require 
further  research  to  complete  analytical  models.  The  current 
modified  exponential  empirical  models  can  be  orders  of 
magnitude  in  error  in  their  density  predictions.  At  higher 
earth  orbits  orders  of  magnitude  error  in  density  is  of 
little  consequence. 3 

The  model  chosen  for  density  and  compositional  effects 
in  the  atmosphere  depends  on  the  importance  of  the  short¬ 
term  cyclic  effects.  Computations  for  orbits  of  months  in 
duration  and  less  than  600  kilometers  in  altitude  will  be 
subject  to  large  errors  using  the  simple  empirical  models 
(Jacchia  1966) .  For  these  short  missions  the  more  accurate 
Jacchia  1977  model  seems  more  appropriate.  (Performance  of 
the  Jacchia  1977  empirical  model  has  been  remarkable  with 
density  errors  rarely  exceeding  two  percent.)  Missions 
lasting  for  decades  may  be  able  to  use  the  computationally 
more  efficient  Jacchia  J60  model  to  predict  very  long  term 
drag  effects. 

The  one  drag  force  dependent  factor  not  yet  discussed  is 
the  satellite-to-atmosphere  relative  velocity.  This 
relative  motion  is  influenced  by  three  gross  factors:  local 


3Above  800  kilometers  the  drag  term  in  the  orbital 
decay  model  used  by  the  Satellite  Control  Facility  in 
Sunnyvale,  California  is  very  small.  Its  main  function  is  a 
catch-all  for  errors  in  other  perturbating  terms. 
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winds,  satellite  tangential  velocity,  and  atmospheric 
rotation. 


Thermospheric  wind  has  not  been  measured.  Mesospheric 
wind  has  been  measured  based  on  the  observation  on  meteor 
trails.  Fleagle  [Ref.  4:p.  82]  reports  meospheric  winds 
exceeding  0.150  kilometers  per  second.  This  may  sound 
large,  but  must  be  compared  with  the  satellite’s  velocity 
(which  is  on  the  order  of  10  kilometers  per  second) . 

Satellite  tangential  velocity  is  a  major  factor  in  the 
satellite-to-atmosphere  relative  velocity.  When  combined 
with  the  rotating  atmosphere  it  accounts  for  the  bulk  of 
relative  motion. 

Most  earth  orbit  calculations  simplify  computation  of 
perturbation- induced  changes  in  the  orbital  elements  by 
assuming  a  fixed  earth.  When  considering  drag  a  non¬ 
rotating  earth  assumption  can  lead  to  a  severe  error.  The 
earth  is  not  fixed.  The  atmosphere  is  pulled  along  with  its 
rotation,  and  lags  the  earth  rotation  by  a  small  fraction. 
Figure  11  shows  the  rotating  atmosphere  effect.  For  low 
inclination  (prograde)  orbits  the  atmosphere  is  actually 
following  the  satellite.  In  high  inclination  (retrograde) 
orbits  the  atmosphere  is  moving  against  the  satellite.  When 
comparing  a  typical  16  revolution  per  day  satellite 
tangential  velocity  in  a  prograde  atmosphere  to  the  same 
satellite  in  a  retrograde  atmosphere,  a  15  percent  variation 
in  relative  velocity  can  occur.  This  is  compared  to  the 


SAME 

ROTATION 


PERPENDICULAR 
ROTATION  CO 


PROGRADE  ORBIT 


POLAR  ORBIT 


OPPOSITE 

ROTATION 


Satellite  to  Atmosphere 


Relative  Rotation 


RETROGRADE  ORBIT 


Figure  11.  Orbital  Inclinations  in  a 
Rotating  Atmosphere 


unknown  (but  probably  less  than  two  percent  effect)  of 
atmospheric  winds. 

A  model  used  by  Parks  [Ref.  l:p.  4)  to  account  for  the 
rotating  atmosphere  separates  the  relative  velocity  of  the 
satellite  in  relation  to  atmosphere  into  two  parts.  The 


model  is  given  by: 
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vrel  ~  vt  1,0  i  J-J-i 

Here  vre^  is  the  satellite  to  atmosphere  relative  velocity, 
v^  is  the  satellite  tangential  velocity,  and  6  is  a  factor 
accounting  for  orbital  orientation  (inclination)  with 
respect  to  a  rotating  atmosphere  (see  Figure  11) . 

Equation  11  in  Parks  [Ref.  l:p.  4]  models  the 
inclination  factor  in  the  following  way: 

r  2 

6  =  [1  -  (^*-)  Awe  cos  i]  (12) 

P 

In  this  model,  6  is  the  rotation  factor,  rp  the  radius  at 
perigee,  Vp  the  satellite  velocity  at  perigee,  Ue  the 
earth's  angular  rotation  rate  (assumed  constant) ,  i  the 
orbital  inclination,  and  A  the  ratio  of  rates  of  rotation  of 
the  atmosphere  to  the  earth's  rotation. 

The  model  for  6  used  in  the  propellant  longevity  model 
neglects  the  effect  of  the  atmospheric  wind,  but  this  should 
introduce  only  a  small  error.  Recent  atmospheric  work, 
Jacchia  [Ref.  7],  suggests  a  significant  gravity  wave 
induced  motion  in  the  upper  atmosphere.  As  more  experimen¬ 
tal  information  on  this  bulk  motion  becomes  available,  the 
incorporation  of  the  gravity  wave  wind  may  become 
increasingly  important. 

The  parameter  a  used  by  the  propellant  longevity  model 
in  the  satellite  velocity  multiplier  submodel  is  probably 
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not  constant.  Atmospheric  theory  suggests  it  may  be 
dependent  on  both  altitude  and  latitude.  The  range  of 
values  for  the  A  parameter  is  assumed  to  be  from  a  maximum 
of  about  1.0  (at  low  altitudes)  to  a  minimum  of  around  0.9 
(at  high  altitudes) . 

The  coding  of  the  theory  in  NSWC  TR  83-243  for  satellite 
fuel  mass  decrement  did  not  compute  6  internally.  It  was 
user  supplied  input  in  the  original  version.  As  part  of  the 
mass  decrement  model  validation  at  the  Naval  Postgraduate 
School,  code  internal  computation  of  6  was  added.  This 
allows  6  to  automatically  adjust  the  satellite-to-atmosphere 
relative  velocity  vrei  for  the  results  predicted  by  the 
theory  as  perturbation- induced  changes  to  the  inclination 
occur.  For  a  long  mission  life  that  exhibits  significant 
inclination  drift,  the  internal  computation  of  should  add 
accuracy  to  the  mass  decrement  submodel.  In  place  of  the 
original  user  input  for  6 ,  the  A  is  now  specified  by  an 
input  file.  The  input  file  default  value  for  A  is  0.9. 

A  model  atmosphere  used  for  drag  computations  must  be 
able  to  correctly  predict  the  effects  on  density,  drag 
coefficient,  and  satellite-to-atmospheric  velocity.  There 
are  many  variations  in  the  atmosphere  that  affect  density. 
The  sum  of  all  these  effects  creates  a  density  bulge  in  the 
atmosphere,  the  location  and  shape  of  which  follow  maximum 
heating. 


The  drag  coefficient  is  atmosphere  composition 
dependent.  Thus  the  model  atmosphere  must  be  able  to 
predict  the  dominant  species  in  the  satellite  environment. 

Satellite-to-atmosphere  relative  velocity  is  a  function 
of  satellite  tangential  speed,  atmospheric  rotation,  and 
orbital  inclination. 

The  Jacchia  J60  model  atmosphere  used  in  the  propellant 
longevity  model  only  includes  the  effects  of  diurnal  heating 
on  the  density  bulge.  There  are  no  computations  of 
atmospheric  composition.  The  newer  Jacchia  J77  model 
incorporates  all  the  known  density  variations  currently 
observed.  The  Jacchia  J77  atmosphere  also  predicts 
thermospheric  composition. 

The  background  concepts  for  the  propellant  longevity 
model  understanding  and  operation  have  been  presented.  The 
next  chapter  describes  the  model  and  its  computer  program. 


This  chapter  describes  the  propellant  longevity  model. 
First,  a  brief  summary  of  the  model  theory  is  presented.  The 
second  part  of  this  chapter  describes  the  coding  of  the 
model  as  installed  at  the  Naval  Postgraduate  School. 
Finally,  the  third  section  gives  the  model  code  compiling 
and  handling  instructions  specific  to  the  system  at  the 
Naval  Postgraduate  School. 

A.  PROPELLANT  LONGEVITY  MODEL  THEORY 

A  detailed  development  of  the  propellant  longevity  model 
is  presented  by  Dr.  Parks  [Ref.  1].  The  propellant 
longevity  model  hinges  on  the  prediction  of  the  fuel  mass 
used  by  a  satellite  to  compensate  for  atmospheric  drag.  The 
mass  decrement  submodel  predicts  the  amount  of  fuel  consumed 
after  each  revolution  of  the  satellite. 

The  mass  decrement  submodel  for  a  low-earth-orbit 
satellite  in  an  oblate  diurnal  atmosphere  is  presented  in 
Reference  1  as  Equation  56.  The  diurnal  density  bulge  is 
the  only  density  and  compositional  atmospheric  variation 
considered  in  the  model. 

The  integrals  of  the  mass  decrement  model  (Equation  56) 
are  next  expressed  in  terms  of  Bessel  functions  of  the  first 
kind  and  order  n  as  defined  by: 
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Here  In  is  the  Bessel  function  of  the  first  kind  and  order 
n,  B  the  inverse  density  scale  height  (assumed  constant) ,  a 
the  semi-major  orbital  axis,  e  the  orbital  eccentricity,  n  = 
0,1,2,  6,  and  E  the  orbital  eccentric  anomaly.  All  of 
these  orbital  mechanics  terms  were  presented  in  Chapter 
II. A.  The  integral  is  to  be  evaluated  over  one  revolution. 

The  Bessel  function  expansion  of  the  mass  decrement 
model  is  given  in  equation  58  of  Reference  1.  This 
particular  expansion  is  central  to  the  model.  Provided 
accurate  values  for  other  multipliers  are  produced  by  the 
model  atmosphere,  (accurate  values  for  atmospheric  density, 
satellite  drag  coefficient,  and  satellite-to-atmosphere 
relative  velocity  are  essential)  greatly  increased  accura¬ 
cies  should  result  in  predicting  the  fuel  consumed  to 
compensate  for  atmospheric  drag.  The  main  draw  back  to  the 
use  of  the  this  type  of  Bessel  function  is  the  known 
slowness  of  its  convergence. 

B.  COMPUTER  CODE  DESCRIPTION 

A  detailed  description  of  the  model  code  and  its  input 
elements  is  contained  in  the  code  and  included  in  Appendix 
A.  An  operational  description  is  presented  here. 

The  computer  code  calculates  the  propellant  mass 
remaining  together  with  the  cumulative  mission  life  after 
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each  revolution  of  the  satellite.  When  propellant  fuel  is 
exhausted,  the  computation  terminates. 

Computer  output  is  available  in  several  forms  and  is 
specified  by  the  user  through  the  use  of  some  control 
parameters.  The  user  can  select  one  of  two  modes;  standard 
mode  or  a  sensitivity  analysis  mode.  The  input  elements  used 
as  control  parameters  for  output  and  mode  control  are: 
maximum  number  of  revolutions  ITER1 ,  sensitivity  analysis 
granularity  PSIZE,  output  print  frequency  ITER2 ,  mode 
control  PSENS,  anomalistic  period  print  control  IET,  number 
of  scheduled  orbit  adjusts  NOA,  scheduled  orbit  adjust  array 
IRV(K),  and  the  amount  of  each  adjust  array  DA ( K )  .  These 
terms  are  discussed  below. 

The  most  basic  form  of  output  available  is  a  tabulation 
of  propellant  remaining  and  accumulated  mission  life  for  20 
specified  satellite-orbital  initial  conditions.  The  20 
initial  conditions  are  specified  in  an  input  file. 
Accumulated  mission  life  is  expressed  in  modified  Julian 
days  (MJD)4  and  revolution  number.  Frequency  of  output 
printing  is  selectable  from  a  maximum  of  once  per  revolution 
to  a  minimum  of  once  per  10000  revolutions  by  user 
specification  or  ITER2 . 

A  cap  on  runaway  computer  time  is  available  by  selection 
of  ITER1.  Satellites  with  small  cross-sectional  areas,  high 

4  A  description  of  the  modified  Julian  day  time 
measurement  is  in  Taff  (Ref.  2:p.  103).  The  modified  Julian 
day  is  reckoned  from  0h  Nov.  17  1858. 


altitudes,  and  large  amounts  of  maneuvering  fuel  can  have 
very  long  orbital  life-times.  By  defining  a  maximum  number 
of  revolutions  to  be  considered  in  the  calculations  through 
the  reasonable  selection  of  ITER1,  unnecessary  expenditure 
of  computer  time  can  be  avoided. 

As  part  of  the  model  verification  a  sensitivity  analysis 
option  was  added  to  the  program.  It  can  be  selected  by 
setting  PSENS  =  1.  The  standard  operating  mode  is  selected 
if  PSENS  =0.  In  the  sensitivity  analysis  mode,  the 
subroutine  SNALYS  creates  an  array  of  elements  selected  by 
the  user  to  be  linearly  scanned  from  a  specified  minimum  to 
maximum  with  a  given  granularity. 

Granularity  is  specified  by  the  selection  of  the  control 
parameter  PSIZE.  A  granularity  of  1/100  (code  specified  as 
PSIZE  *  100)  will  cause  100  executions  of  the  propellant 
longevity  code  to  complete  fuel  exhaustion.  This 
specification  can  lead  to  excessive  computer  execution  times 
for  high  altitude  satellites.  A  granularity  of  1/5  will 
cause  five  executions  of  the  code  to  fuel  exhaustion. 
However,  only  five  iterations  from  a  given  minimum  to  a 
given  maximum  may  be  too  coarse  to  reveal  the  true  trend  of 
the  satellite  mission  life. 

Each  time  the  code  is  sequentially  executed  in  the 
sensitivity  analysis  mode,  the  elements  selected  by  the  user 
for  analysis  are  incremented  linearly.  The  output  in  the 
sensitivity  analysis  mode  is  a  table  of  mission  lifetimes, 


resulting  from  the  initial  values  of  input  elements  not 
selected  for  analysis,  and  the  current  value  of  the  elements 
selected  for  analysis. 

The  output  always  echoes  the  initial  input  element 
conditions  as  read  from  a  user  supplied  input  file.  Further 
choice  of  output  and  mode  control  parameters  allow  the  user 
to  tailor  the  output. 

Three  more  input  control  elements  allow  for  the 
specification  by  revolution  number  of  up  to  ten  adjusts  to 
the  semi-major  axis.  The  element  IRV(K)  is  an  array  whose 
elements  specify  the  orbital  revolution  number  for  axis 
adjust.  The  element  DA(K)  specifies  the  amount  of  change  in 
the  semi-major  axis  in  kilometers  for  each  scheduled 
adjustment.  Finally,  NOA  specifies  the  total  number  of 
orbital  adjusts  to  be  performed.  If  the  user  does  not 
desire  any  changes  to  the  semi-major  axis,  NOA  must  be  set 
to  0  and  no  values  for  IRV(K)  and  DA(K)  input  to  the  file 
stream. 

An  additional  twenty  input  elements  must  be  specified  by 
the  user  in  the  input  file.  Five  of  these  are  satellite 
physical  parameters,  five  are  atmosphere  specific 
parameters,  and  ten  are  the  Brouwer  mean  orbital  element  set 
of  the  NAVSPASUR  system. 

The  five  satellite  physical  parameters  are;  propulsion 
motor  specific  impulse  (SPI  in  seconds) ,  initial  mass  of  the 
satellite  including  fuel  (WT  in  kilograms),  initial  mass  of 
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the  maneuvering  fuel  (W  in  kilograms) ,  satellite  atmospheric 
drag  coefficient  (CD) ,  and  satellite  cross-sectional  area  (A 
in  square  kilometers) . 

Motor  specific  impulse  SPI  can  range  from  220  seconds 
for  monopropellants  to  3000  seconds  for  ion  drives.  The 
default  value  is  230. 

The  initial  satellite  mass  WT  depends  on  the  satellite 
mission  and  design.  A  typical  communications  satellite  has 
a  mass  of  about  3500  kilograms.  An  SDI  satellite  may  be 
considerably  more  massive.  The  default  value  for  the 
satellite  mass  is  9100  kilograms. 

Initial  maneuvering  fuel  mass  W  of  the  satellite  depends 
on  the  design  life  of  the  satellite  as  well  as  it's 
ballistic  coefficient.  Low  altitude  satellites  need  at 
least  one  percent  of  their  initial  mass  as  maneuvering  fuel 
if  mission  lives  of  years  are  to  be  achieved.  The  default 
value  for  initial  maneuvering  fuel  is  100  kilograms. 

The  drag  coefficient  CD  is  a  measure  of  the  slowing 
power  of  the  atmosphere  on  a  given  satellite.  In  the 
thermosphere  it  is  very  atmosphere  composition  dependent.  A 
value  of  2.2  has  been  settled  on  by  orbital  analysts  for 
altitudes  of  200  to  400  kilometers.  A  value  of  at  least  3.0 
must  be  assigned  at  altitudes  above  500  kilometers, 
(corresponding  to  the  helium  atmosphere) .  The  default  value 
for  the  drag  coefficient  is  2.2. 


Satellite  cross-sectional  area  is  again  very  design 
dependent.  The  echo  series  satellites  had  cross-sections 
measuring  in  hundred's  of  square  meters.  Some  SDI 
satellites  may  also  have  very  large  cross-sections.  The 
default  for  cross-sectional  area  is  50.0*10“6  square 
kilometers,  corresponding  to  50  square  meters. 

The  five  atmosphere  specific  parameters  are:  atmosphere- 
to-earth  relative  rotation  rate  D,  10.7  cm  solar  flux  F107 
described  in  Chapter  II,  ninety-day  average  solar  flux  FBAR, 
geomagnetic  activity  index  AKP,  and  the  time  of  vernal 
equinox  passage  TVE. 

The  atmosphere-to-earth  relative  rotation  rate  D  is  the 
factor  explained  at  the  end  of  Chapter  II.  The  default 
value  is  0.9. 

The  decimetric  solar  flux  and  the  average  decimetric 
solar  flux  F107  and  FBAR  are  again  explained  in  Chapter  II. 
Default  values  of  100  are  assigned  to  each. 

The  geomagnetic  index  AKP  is  discussed  in  Chapter  II. 
It  ranges  from  0  to  7  and  is  set  at  a  default  value  of  2.00. 

With  the  current  atmospheric  model  Jacchia  J60,  both 
FBAR  and  AKP  are  read  as  input  but  not  used  for  the  density 
computation.  As  seen  from  a  sensitivity  analysis,  these 
parameters  currently  have  no  effect  on  mission  life.  The 
Jacchia  J77  model  atmosphere  uses  these  now  stranded  input 
parameters.  The  future  incorporation  of  the  Jacchia  J77 
atmosphere  model  will  be  able  to  use  these  input  elements. 


The  ten  NAVSPASUR  orbital  elements  are  the  final  input 
elements  in  the  standard  input  file.  They  are  all  explained 
in  the  orbital  mechanics  section,  and  in  the  code  included 
in  Appendix  A.  Briefly,  they  are: 

UJD:  time  of  epoch  in  modified  julian  days 

ETU:  initial  mean  anomaly 

GO:  argument  of  perigee 

HO:  right  ascension  of  the  ascending  node 

ES:  eccentricity 

XI:  inclination 

AO:  semi-major  axis 

RND:  rate  of  change  of  mean  motion 

ESD:  eccentricity  decay  rate 

ADOT:  semi-major  axis  decay  rate. 

The  first  seven  of  these  are  orbital  elements  and  the  final 
three  are  the  results  of  perturbations.  The  NAVSPASUR 
orbital  element  set  is  the  result  of  orbital  measurements  of 
actual  satellites. 

The  final  block  of  input  in  the  input  file  is  the 
sensitivity  analysis  specification  array.  It  must  be 
included  only  when  the  sensitivity  analysis  mode  is 
specified.  The  array  consists  of  20  lines  used  to  make  up 
the  array  of  elements  to  be  scanned  by  the  program  when  the 
sensitivity  analysis  mode  is  selected. 

The  sensitivity  analysis  specification  array  is  divided 
into  four  columns.  The  first  column  contains  the  input 
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element  variable  name. 


These  are  the  same  names  as 


specified  in  the  satellite  physical,  atmosphere  specific, 
and  orbital  elements  input  discussions. 

The  second  column  is  a  control  parameter.  When  set  to 
1.0,  the  corresponding  input  element  is  included  in  the 
sensitivity  analysis.  Columns  three  and  four  are  in  a 
Fortran  F16.9  format  and  they  specify  the  minimum  and 
maximum  range  of  the  element  to  be  scanned.  Defaults  for 
the  minimum  and  maximum  have  been  provided  in  specific  cases 
corresponding  to  logical  and  physical  ranges  of  the  input 
elements.  In  some  cases  zeros  appear  in  the  minimum  and 
maximum  columns.  These  zeros  correspond  to  input  elements 
not  physically  constrained.  An  example  of  an  input  element 
not  physically  constrained  is  the  time  of  epoch  (UJD) . 

If  a  1.0  is  entered  in  column  two  of  a  particular  row  of 
the  sensitivity  analysis  specification  array,  that  input 
element  is  included  in  the  analysis.  Any  number  of  the 
input  elements  may  be  selected  for  simultaneous  analysis. 
When  selecting  a  multivariable  sensitivity  analysis,  the 
user  must  understand  the  interactions  of  the  incrementing 
variables  and  ensure  that  a  solution  exists  in  the  mission 
life  range  being  considered. 

When  a  variable  is  included  in  the  sensitivity  analysis, 
the  incrementing  values  override  the  initial  conditions. 
All  input  elements  not  user  selected  for  analysis  (those 
with  zeros  in  the  second  column)  are  reset  to  their  initial 
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values  after  each  iteration  during  the  analysis.  Appendix  B 
contains  sample  input  element  files  and  their  corresponding 
outputs . 

Input  elements  and  output  forms  have  been  explained. 
The  program  flow  is  discussed  next.  The  16  subroutines  and 
eight  functions  comprising  the  propellant  longevity  model 
are  explained  in  Appendix  A  in  the  Fortran  code.  The 
routine  PLEP  initializes  the  system,  including  inputting 
program  control  parameters  and  initial  values  for  all 
elements.  If  the  sensitivity  analysis  mode  is  selected, 
subroutine  SNALYS  generates  the  array  to  be  scanned 
linearly,  complete  with  the  necessary  increments  for  each 
iteration.  Analysis  mode  causes  completion  of  all 
computations  to  propellant  exhaustion  PSIZE  times.  The 
routines  SNALYS  and  PLEP  are  completed  only  once  even  in 
analysis  mode. 

Program  control  is  transferred  from  the  initialization 
routine  PLEP  to  the  subroutine  KOZAI,  the  main  driver  and 
bookkeeper  of  the  system.  The  routine  KOZAI  first  causes 
the  NAVSPASUR  Brouwer  mean  element  set  to  be  converted  into 
a  Keplerian  (Kozai-like)  mean  orbital  set,  by  calling  the 
subroutines  BRAUER,  PERIOD,  and  MEAN. 

Subroutine  MEAN  calls  functions  XSPA-XSPM  to  convert  the 
Brouwer  osculating  elements  generated  in  BRAUER  to  Keplerian 
orbital  elements.  The  Keplerian  elements  are  then  returned 
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to  the  subroutine  KOZAI.  In  standard  mode  the  mean  element 
set  is  given  as  output. 

All  of  the  above  computations  have  simply  conditioned 
the  input  elements  to  conform  to  the  requirements  of  the 
propellant  longevity  model.  The  routine  KOZAI  now  causes 
calculation  of  the  fuel  used  to  compensate  for  atmospheric 
drag  after  each  revolution.  The  computation  is  continued 
until  the  fuel  is  exhausted  or  the  maximum  number  of 
revolutions  (ITER1)  is  reached.  When  the  computation 
terminates,  program  control  is  transferred  back  to  PLEP. 

To  calculate  the  drag-compensating  fuel  required  for 
each  revolution,  KOZAI  calls  the  subroutines  PERIOD,  GEOP, 
THRST,  and  DRAG. 

The  subroutine  GEOP  calculates  the  change  in  orbital 
elements  due  to  the  aspherical  earth.  It  uses  the  J2,  J3, 
and  J4  terms  in  a  Legendre  polynomial  expansion,  as 
explained  in  Chapter  II.  Aspherical  earth  calculations  are 
performed  once  each  revolution. 

The  subroutine  THRST  checks  each  revolution  for  a 
scheduled  orbital  adjust.  If  a  change  to  the  semi-major 
axis  is  scheduled,  it  calculates  the  change  to  the  axis,  the 
argument  of  perigee,  the  eccentricity,  and  the  mean  anomaly. 
It  also  computes  the  amount  of  maneuvering  fuel  required  to 
perform  the  adjustment.  All  changes  are  returned  to  the 
subroutine  KOZAI . 


The  subroutine  DRAG  calculates  the  fuel  mass  needed  to 
overcome  the  atmospheric  drag.  This  calculation  is 
accomplished  by  calling  four  subroutines:  SOLOC,  SATLOC, 
SALT,  and  JAC60 . 

The  routine  SOLOC  computes  the  right  ascension  and 
declination  of  the  sun  for  the  model  atmosphere  diurnal 
density  bulge.  The  routine  SATLOC  computes  the  right 
ascension,  declination  and  geomagnetic  latitude  of  the 
satellite.  This  is  accomplished  by  first  calling  the 
subroutine  POSVEL  to  calculate  satellite  position  and 
velocity.  The  routine  POSVEL  in  turn  calls  the  subroutine 
NWTRPH  to  solve  Kepler's  equation  (see  Chapter  II).  The 
routine  SALT  computes  the  satellite  altitude. 

Atmospheric  density  is  computed  by  the  subroutine  JAC60 
based  on  the  Jacchia  J60  model  atmosphere  discussed  in 
Chapter  II.  A  maximum,  minimum,  and  mean  density  that  the 
satellite  will  encounter  in  the  oblate  diurnal  atmosphere 
are  calculated. 

The  mass  decrement  equation  (Equation  58  of  Reference  1) 
is  contained  in  the  subroutine  DRAG.  The  function  MMBSIR  is 
used  to  perform  the  Bessel  function  expansion.  The  routine 
DRAG  returns  the  fuel  mass  decrement  due  to  atmospheric  drag 
computation  to  KOZAI.  The  routine  KOZAI  then  decrements  the 
remaining  fuel  by  this  amount  and  checks  for  fuel  exhaustion 
and  maximum  orbital  number.  If  fuel  is  less  then  zero,  or 
the  current  revolution  number  equals  ITER1 ,  calculations  are 


terminated  and  control  is  returned  to  PLEP.  In  standard 
mode,  program  execution  now  halts. 

In  sensitivity  analysis  mode,  all  input  elements  not 
selected  for  analysis  are  returned  to  their  initial  values, 
elements  selected  for  analysis  are  linearly  incremented,  and 
K02AI  is  again  called  by  PLEP.  The  subroutine  KOZAI  then  is 
executed,  as  previously  described,  with  the  new  values  for 
the  selected  elements.  Upon  fuel  exhaustion  or  reaching  the 
maximum  revolution  number,  KOZAI  terminates  calculation  and 
passes  program  control  back  to  PLEP.  This  process  continues 
until  the  elements  user  selected  for  analysis  have  been 
incremented  from  their  minimum  to  maximum.  When  the  maximum 
value  of  the  elements  user  selected  for  analysis  is  reached, 
program  execution  terminates  in  the  sensitivity  analysis 
mode. 

The  tabular  output  in  the  analysis  mode  is  printed  line- 
by-line  in  KOZAI  upon  reaching  fuel  exhaustion  or  maximum 
revolution  number.  The  result  is  a  table  of  data,  PSIZE 
long,  that  represents  the  linear  scan  of  the  specified 
elements.  The  first  line  is  the  mission  life  corresponding 
to  the  minimum  value  of  the  scanned  elements.  Each 
subseguent  line  represents  the  total  mission  life  resulting 
from  an  increment  of  the  elements  selected  for  analysis  and 
the  initial  conditions  of  the  elements  not  selected.  The 
final  line  is  the  mission  life  resulting  from  selected 


elements  at  their  maximum  value  and  non-selected  elements  at 
their  initial  conditions. 

Mission  life  is  expressed  in  revolution  number  and 
modified  Julian  days.  The  last  column  of  the  tabular  output 
in  analysis  mode  represents  the  current  value  of  the 
incremented  input  element.  When  more  than  one  element  is 
analyzed,  the  last  column  is  the  current  fraction  of  scan,  a 
real  number  between  0  and  1. 

The  tabular  output  can  easily  be  converted  into  an  input 
file  for  a  graphics  routine  by  deleting  the  first  98  lines 
of  output. 

C.  MODEL  CODE  OPERATION 

This  section  describes  the  operation  of  the  computer 
code  at  the  Naval  Postgraduate  School  on  the  IBM  3  03  3 
computer  under  VMS  using  the  VS  FORTRAN  compiler. 

The  propellant  longevity  model  must  be  compiled  in  the 
VM  mode  using  the  automatic  precision  increase  option.  If 
the  propellant  longevity  code  is  in  filename  SDI  and 
filetype  FORTRAN  (the  filetype  must  be  FORTRAN) ,  the  correct 
form  for  the  compilation  command  is: 

FORTVS  SDI  (AUTODBL(DBLPAD4) ) 

The  listing  generated  by  this  compilation  takes  up  about  20 
percent  of  the  user’s  "A"  disk. 

Prior  to  program  execution,  the  input  and  output  files 
must  be  specified.  If  the  input  file  is  in  filename  PLEP01 
and  filetype  DATA2,  then  the  correct  form  of  the  command  is: 
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FILEDEF  5  DISK  PLEP01  DATA2  A1 

This  informs  the  computer  where  to  look  for  input  during 
execution  of  the  model  code.  If  hard  copy  of  the  output  is 
desired,  it  can  be  written  to  a  file  instead  of  the  screen 
by  using  the  FILEDEF  command.  If  the  file  name  is  SDI  and 
the  desired  type  is  DATA,  the  correct  command  is: 

FILEDEF  6  DISK  SDI  DATA  A1 

After  execution  of  the  propellant  longevity  code,  the 
results  will  be  in  SDI  DATA  on  the  user's  "A"  disk. 

To  execute  the  model  listing  generated  upon  compilation, 
the  command  RUN  SDI  is  typed  (this  assumes  that  SDI  is  the 
filename  of  the  propellant  longevity  program) . 

This  program  can  have  very  long  execution  times.  A 
general  rule  of  thumb  is  that  execution  time  takes  less  than 
PSIZE  *  ITERl/10 , 000  minutes.  Roughly,  the  computation 
takes  100  microseconds  for  each  revolution.  For  example,  if 
PSIZE  is  100  and  ITER1  is  10,000  (corresponding  to  a  five 
year  lifetime)  then  the  IBM  3033  will  take  200  minutes  to 
complete  the  calculations.  The  user  must  carefully  choose 
small  values  of  PSIZE  and  ITER1  on  the  initial  investigation 
of  elements. 

The  simplest  graphical  routine  to  construct  plots  of 
tabular  data  is  EASYPLOT.  The  EASYPLOT  routine  was  used  to 
generate  all  the  plots  shown  in  Chapter  IV. 

In  summary,  the  steps  of  program  operation  are: 

(1)  Compile  the  program  using  the  automatic 
precision  increase  option. 
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(2)  Select  the  mode  of  operation  and  input  element 
values  by  manipulation  of  the  input  file. 

(3)  Define  input  and  output  files. 

(4)  Execute  the  program  by  typing  the  file  name. 

The  next  chapter  presents  the  results  of  the  input 
element  sensitivity  analysis  performed  to  test  the  model 
reasonableness  during  verification. 


IV. 


This  chapter  presents  the  results  of  the  propellant 
longevity  model  verification  at  the  Naval  Postgraduate 
School.  Standard  model  verification  tests  the  model  against 
three  questions.  First,  does  the  model  address  the  problem? 
Second,  are  the  model  predictions  reasonable?  Third,  how  do 
the  model  predictions  compare  to  real-world  observations? 

This  verification  addresses  only  the  first  two  criteria. 
Observational  data  for  comparison  with  the  model  predictions 
is  not  available  in  an  unclassified  state.  The  quantitative 
analysis  of  model  predictions  is  thus  not  addressed  in  this 
verification. 

The  model  does  address  the  initial  problem.  It  predicts 
the  fuel  mass  decrement  of  a  low-earth-orbit  satellite  doing 
intrack  micro-thrusting  to  overcome  atmospheric  drag. 

The  reasonableness  test  is  the  subject  of  the  remainder 
of  this  chapter.  An  effective  way  to  test  the  model 
performance  is  to  examine  model  predictions  resulting  from 
the  full  "intended  use"  range  of  the  input  variables. 
Furthermore,  if  a  plot  is  made  of  predicted  results  as  a 
selected  input  element  is  incremented  through  its  range, 
prediction  trends  can  be  found  and  examined.  This  procedure 
provides  a  qualitative  analysis  of  model  sensitivity  to  the 
selected  input  element. 
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To  aid  in  these  tests  and  trend  studies,  a  sensitivity 
analysis  feature  was  incorporated  into  the  program  code  of 
the  model,  as  previously  described.  Graphical  results  of 
each  single  variable  sensitivity  analysis  are  presented  for 
the  input  elements. 

The  format  of  the  following  figures  are  all  related. 
They  represent  plots  of  satellite  mission  life  in  modified 
Julian  days  (here  after  referred  to  simply  as  "days")  versus 
a  specified  input  element.  All  figures  have  the  Y-axis  as 
mission  life  and  the  X-axis  as  the  variable  selected  for 
analysis.  In  each  case  the  range  of  values  for  the  input 
element  is  selected  to  remain  within  reasonable  values  of 
the  input  variable,  as  discussed  in  Chapter  III.  The 
granularity  in  most  cases  is  1/25  (PSIZE  =  25)  .  In  some 
cases  it  is  necessary  to  decrease  the  granularity  to  l/ioo 
(PSIZE  *  100)  to  see  the  actual  trends. 

The  default  values  discussed  in  Chapter  III  are  used  for 
all  input  elements  not  selected  for  analysis  (except  where 
noted) .  They  correspond  to  a  near-polar,  retrograde  orbit 
of  an  actual  satellite  supplied  as  sample  input  by  NSWC. 
The  model  predicted  propellant  longevity  using  all  default 
conditions  and  10  kilograms  of  fuel  is  13.74  days.  When  the 
default  values  and  100  kilograms  of  fuel  are  used,  the 
resulting  propellant  longevity  is  150.21  days.  An  *  symbol 
on  the  plotted  curves  indicates  the  default  value  of  the 


i 

$ 

r 

t 

r 

i 

k 

k 

k 

k 

K 

a 


analyzed  input  element.  The  default  value  names  are  again 


summarized  here.  The  units  can  be 

Satellite  Dhvsical  Darameters: 

found  in  Chapter 

SPI  (motor  specific  impulse) 

230 

WT  (satellite  initial  mass) 

9100 

W  (initial  fuel  mass) 

100  and  10 

CD  (drag  coefficient) 

2 . 2  and  3 . 0 

A  (cross-sectional  area) 

. 00005 

AtinasBfaeris.  ..sKssifis  car  a, meters,;.. 

D  (atm.  rotation  factor) 

F107  (decimetric  solar  flux) 

FBAR  (average  F107) 

AKP  (geomagnetic  index) 

TVE  (time  of  vernal  equinox) 

FAYS PAS UR  qrbital  elements; 

UJD  (time  of  epoch) 

ETU  (mean  anomaly) 

HO  (R. A.  of  ascending  node) 

GO  (argument  of  perigee) 

ES  (eccentricity) 

XI  (inclination) 

RND  (time  rate  of  mean  motion)  .000783912 
ESD  (eccentricity  decay  rate)  -.239590000 
A0  (semi-major  axis)  1.06 

ADOT  (semi-major  axis  decay)  -701.8652 


1.0 

100 

100 

2.0 

43222.7382 

44619.98716775 
150.0000 
95.0000 
200.0000 
. 0100000 
95.0000 
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The  input  element  sensitivity  analyses  are  divided  into 
related  sets.  All  the  sets  are  discussed  in  the  same  basic 
format.  Each  analysis  is  single-variable  and  arranged  as 
follows: 

(1)  The  expected  trend  of  real-world  mission  life 
resulting  from  the  selected  variable 
incrementing  through  its  range. 

(2)  The  expected  mission  life  trend  predicted  by 
the  model  as  a  result  of  the  selected 
variable  incrementing  through  its  range. 

(3)  The  actual  mission  life  results  of  the 
selected  variable  incrementing  through  its 
range . 

(4)  A  discussion  of  the  selected  variable 
sensitivity  analysis. 

The  explanations  for  the  plotted  trends  are  based  on 
expected  behavior  in  agreement  with  current  theory.  An 
understanding  of  the  actual  mechanisms  producing  the 
observed  trends  will, be  improved  with  further  research. 

For  a  given  set  of  satellite-orbital  configurations,  the 
satellite  mission  life  will  be  shortest  for  satellites 
traveling  through  the  most  dense  atmosphere.  Thus  for  a 
given  input  variable  sensitivity  analysis,  the  range  of 
values  corresponding  to  satellite  travel  through  the  most 
dense  atmospheric  conditions  will  produce  the  shortest 
mission  life.  The  range  of  values  of  the  same  input 
variable  corresponding  to  satellite  travel  through  the  least 
dense  atmosphere  will  likewise  produce  the  longest  mission 
life.  The  orbital-perigee  atmospheric-density-bulge 


encounter  geometry  is  thus  critical  in  mission  life 
analysis. 

The  initial  sensitivity  plot  of  an  input  variable  is 
accomplished  with  10  kilograms  of  fuel.  In  many  cases,  true 
geopotential-induced  orbital-drift  trends  are  not  seen  with 
small  amounts  of  drag-compensating  fuel.  To  verify  trends 
caused  by  geopotential-induced  orbital  drift,  the  initial 
compensating  fuel  is  increased  to  100  kilograms.  The 
differences  in  the  10  kilogram  and  100  kilogram  mission  life 
plot  trends  can  then  be  attributed  to  geopotential  effects 
and  the  movement  of  the  center  of  most  direct  solar  heating 
during  the  mission. 

The  133.14  day  nominal  mission  life  resulting  from  the 
default  values  of  the  input  elements  and  100  kilograms  of 
drag  compensating  fuel  is  almost  a  half-year.  During  this 
half-year,  the  earth  is  revolving  around  the  sun  causing  the 
noon-sun  atmospheric  position  to  increase  in  longitude  by 
about  one  degree  per  day.  This  greatly  complicates  trend 
analysis  when  using  the  100  kilogram  fuel  case  alone.  It  is 
thus  necessary  to  include  the  10  kilogram  fuel  case.  The 
13.74  nominal  mission  life  of  the  10  kilogram  fuel  case  is 
short  enough  so  that  it  shows  little  longitudinal  movement 
of  the  most-direct-sun  atmospheric  position. 

In  many  of  the  sensitivity  analysis  plots,  the  predicted 
mission  life  curve  has  a  "knee.”  This  "knee"  corresponds 
to  the  range  of  the  selected  input  variable  that  produces 


changing  behavior  in  the  mission  life  curve.  The  location 
of  these  ranges  has  significance  to  design  research. 

The  input  variables  are  divided  into  four  categories 
based  upon  their  effect  on  mission  life.  The  first  category 
contains  the  input  variables  that  affect  the  atmospheric- 
bulge  orbital-perigee  encounter  geometry.  The  second 
category  contains  the  input  variables  that  affect  the 
slowing  effectiveness  of  the  atmospheric  density  bulge  on 
the  particular  satellite  configuration.  The  third  category 
contains  the  input  variables  that  affect  the  depth  of 
penetration  of  the  satellite  into  the  atmosphere  at  perigee. 
Finally,  the  fourth  category  contains  the  input  variables 
that  have  no  affect  at  all  on  mission  life.  As  will  be 
seen,  predicted  mission  life  is  most  sensitive  to  the  third 
category  (factors  affecting  atmospheric  penetration  at 
perigee) . 

Factors  affecting  the  atmospheric-density-bulge  orbital  - 
perigee  encounter  geometry  are:  the  time  of  vernal  equinox 

passage  TVE,  the  time  of  epoch  UJD,  the  argument  of  perigee 
GO,  the  right  ascension  of  the  ascending  node  ho,  and  the 
inclination  XI.  The  time  of  vernal  equinox  provi  ies  i  real 
time  reference  for  earth  in  its  orbit  relative  *  •'he  ;mi 
When  combined  with  the  time  of  epoch  .sometime-;  a,  lei  *  he 
t  i  me  o  f  pe  r  i  gee  passage  )  ,  these  •wo  f  a-  t  ■  -  r  s  Iff  m  •  '  *■ 

position  of  the  noon  sun  in  right  as-erisc  r  m  l  le  "»•  • 

relative  to  the  atmosphere.  I  he  d»-  f  a  u  1  •  i  ;  c  ■  •  i  *  ’  * 


model  verification  position  the  initial  noon  sun  at 
approximately  300  degrees  longitude  and  20  degrees  south 
latitude  relative  to  the  celestial  meridian  containing  the 
first  point  of  aries.  The  right  ascension  of  the  ascending 
node  and  the  argument  of  perigee  combine  to  position  the 
point  of  perigee  in  longitude  relative  to  the  first  point  of 
aries.  The  orbital  inclination  positions  the  point  of 

perigee  in  latitude  (or  declination) .  Actual  positioning  of 
the  satellite  in  the  above  orbit  is  accomplished  by 
specifying  the  mean  anomaly.  Ail  these  orbital  mechanics 
terms  have  been  explained  in  Chapter  I  I. A. 

The  expected  real-world  behavior  of  satellite  mission 
life  as  a  result  of  varying  the  time  of  vernal  equinox  or 
time  of  epoch  is  discussed  next.  The  time  of  vernal  equinox 
TVE  serves  to  set  a  real-time  reference  for  posit loninq  ot 
the  earth  in  its  orbit  relative  to  the  sun.  When  the  time 
of  epoch  UJD  is  specified,  the  position  of  the  point  of  rh<* 
most  direct  atmospheric  heatinq  i noon  sum  in  the  atmosphere 
is  set.  The  atmospheric  den  sit  /  -bulge  follows  the  1  <  ic  a  M  ■  -  n 
of  the  most  direct  solar  heating.  The  result  of  P.’F  and  ' 
varying  in  relation  ♦  o  each  of  her  is  a  sea  son  a  1  posit  i  •  m  i  rvj 
of  t^e  earth  i  ri  ;  ♦  s  >t  h  1 1  ar-urvl  t  fie  sun  at  ♦  he  *  •  me  ? 
initial  perigee  pass  age  !•*.**  •  f  he  earth's  r  -*t  :  •  -  na  ;  axis 

t  i  i  t  ,  t  fie  point  ■  f  most  1  ,  r  e  t  ;  a  r  he  a»  1  n  j  m  • 

1  onq  i  t  ude  aril  1  a  *  .  ♦  j  le  i  »  r.e  earth  m.  e-.  ;  ,  ♦  ■ .  ?t  .  * 


around  the  sun.  Figures  3  and  4  in  Chapter  II. A  aid  this 
visualization . 

The  earth  and  solar  cycles  effecting  the  atmospheric 
density-bulge  combine  to  define  its  size  and  shape.  These 
cycles  are  discussed  in  Chapter  II. B,  and  stated  again  here. 
They  are:  diurnal,  semi-annual,  seasonal-latitudinal,  solar 
rotational,  and  solar  activity  related.  The  Jacchia  J60 
model  atmosphere,  as  used  in  the  propellant  longevity  model, 
incorporates  only  the  diurnal  cycle. 

Once  the  position,  size,  and  shape  of  the  atmosphere  is 
set,  the  mission  life  of  a  drag  compensating  satellite 
orbiting  within  it  can  be  developed.  For  such  a  satellite 
in  a  specific  initial  orbit,  the  following  real-world  trend 
in  mission  life  (propellant  longevity)  should  be  observed. 
The  trend  should  be  cyclic,  resulting  from  the  superposition 
of  all  the  atmospheric  heating  cycles.  Mission  lives  will 
be  shortest  for  orbits  having  their  perigee  in  the  most 
dense  atmospheric  bulges.  The  longer  mission  lives  will  be 
observed  for  satellites  having  their  orbits  farthest  away 
from  the  density-bulge.  Thus  the  plots  of  satellite  mission 
life  will  show  peaks  and  valleys  consistent  with  the  above 
d 1 scuss l on . 

In  the  propellant  longevity  model,  the  diurnal  cycle  is 
the  only  atmospheric  heating  cycle  that  is  incorporated. 
Positioning  of  the  earth  in  its  orbit  is  accomplished  by  the 
combination  of  the  time  of  vernal  equinox  and  * i me  of  epoch 
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in  the  subroutine  SOLOC.  As  discussed  in  Chapter  III,  the 
subroutine  SOLOC  positions  the  sun  by  right  ascension  and 
declination  over  the  earth's  atmosphere.  The  routine  DRAG 
considers  the  diurnal  density  bulge  to  be  centered  under  the 
solar  position.  The  density  profile  of  the  atmospheric 
density-bulge  is  generated  in  the  routine  JAC60  based  on  the 
value  of  the  input  element  F107 . 

The  propellant  longevity  model  considers  the  atmosphere 
to  be  oblate,  diurnal,  and  revolving  around  the  sun  with  a 
23.5  degree  rotational-axis  tilt.  The  expected  model- 
predicted-mission-life  trend  as  is  cyclic  as  TVE  and  UJD  are 
varied.  Every  year  should  be  an  exact  replica  of  every 
other  year.  Two  peaks  and  two  valleys  in  the  yearly 
mission-life  plot  should  be  observed.  The  deepest  valley 
occurs  when  the  point  of  perigee  is  most  closely  centered  in 
the  atmospheric  density  bulge.  With  a  low  eccentricity  (e  * 
.01)  orbit  considered,  a  much  shallower  valley  will  occur 
about  six  months  later.  This  shallower  valley  corresponds 
to  the  point  of  apogee  being  most  closely  centered  in  the 
atmospheric  density-bulge. 

Two  peaks  in  plotted  mission  life  should  occur 
corresponding  to  the  orbital  pane  being  most  perpendicular 
to  the  earth-sun  line.  The  peaks  should  be  six  months 
apart,  of  approximately  equal  magnitude,  and  interspaced 
between  the  mission-life  valleys.  As  the  eccentricity  of 
the  orbit  is  increased,  the  apogee-associated  valley  will 
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become  shallower  while  the  perigee-associated  valley  becomes 
deeper. 

Figures  12,  13,  14,  and  15  present  the  plots  of  the 
sensitivity  analysis  of  TVE.  In  Figure  12,  TVE  has  a  range 
of  1000  days,  from  38422  to  39422  days.  The  initial  fuel  is 
10  kg.  The  resultant  predicted  mission  life  ranges 
cyclically  between  11  and  18  days. 

Figure  13  has  an  incremented  range  for  TVE  of  38422  to 
38522,  or  100  days.  Initial  fuel  is  10  kg.  The  resulting 
predicted  mission  life  again  ranges  cyclically  between  11 
and  18  days. 

In  Figure  14  the  width  of  the  incremented  range  is 
raised  to  10,000  days  (about  27  years).  The  minimum  value 
of  TVE  is  38422  and  the  maximum  is  49422.  Again  10  kg  of 
fuel  is  used.  The  results  are  cyclic  between  11  and  18  days 
until  TVE  becomes  greater  than  time  of  epoch  UJD.  The 
results  at  the  high  end  of  TVE  in  this  plot  are  thus 
unreasonable.  The  time  of  vernal  equinox  passage  chosen 
must  occur  before  the  satellite  is  launched  requiring  TVE  to 
be  less  then  UJD. 

The  TVE  is  again  incremented  through  a  small  range  of  20 
days  in  Figure  15.  Initial  fuel  is  100  kg.  The  resulting 
range  of  mission  life  predicted  is  from  150  to  153  days. 
Note  the  absence  of  the  cyclic  trend. 

Figures  16  and  17  present  the  results  of  the  sensitivity 
analysis  plots  of  the  time  of  epoch  UJD  with  a  range  of  1000 
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Figure  12.  Mission  Life  versus  Time  of  Vernal 
Equinox  (1000  Days) 
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Figure  13.  Mission  Life  versus  Time  of 
Vernal  Equinox  (100  Days) 
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Figure  17. 


Mission  Life  versus  Time  of  Epoch 
(100  Days  and  100  kg  of  Fuel) 


as 


days.  Note  the  similarity  to  the  TVE  plots, 

theoretically  expected.  Figure  16  shows  the  10  kg  fuel 
case,  and  Figure  17  the  100  kg  fuel  case. 

The  trends  observed  in  Figures  12  through  17  are  in 
agreement  with  expected  model  predictions.  The  plotted 
trends  may  have  significant  differences  from  real-world 
trends  due  to  the  limitations  of  the  Jacchia  J60  model 
atmosphere.  It  should  be  noted  that  the  default  values  of 
TVE  and  UJD  produce  a  s.olar  location  approximately  3  0 
degrees  longitude  past  the  winter  solstice  in  the  southern 
hemisphere  heading  toward  the  spring  (vernal)  equinox.  The 
approximate  latitude  of  the  default  solar  position  is  20 
degrees  south. 

The  right  ascension  of  the  ascending  node  HO  positions 
the  orbital  plane  with  reference  to  the  celestial  meridian 
containing  the  first  point  of  aries.  When  the  argument  of 
perigee  GO  and  the  inclination  XI  are  specified,  the  point 
of  perigee  is  positioned  in  longitude  and  latitude.  This 
geometry  is  shown  on  Figure  4  in  Chapter  II. A. 

The  real-world  trend  in  mission  life  of  a  drag 
compensating  satellite,  as  a  result  of  varying  HO,  GO,  and 
XI  is  cyclic.  The  cycle  repeats  every  360  degrees  for  HO 
and  GO,  and  every  180  degrees  for  XI.  The  observed  pattern 
of  peaks  and  valleys  will  be  highly  variable,  depending  on 
the  values  of  these  three  factors. 
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Model  predicted  mission  life  trends  due  to  varying  HO, 
GO,  and  XI  through  their  ranges  should  closely  parallel 
real-world  expectations.  This  close  parallel  between  model 
predictions  and  observation  assumes  the  atmospheric  density- 
bulge  has  been  properly  positioned  and  described. 

The  default  values  of  95  and  200  degrees  for  the 
inclination  and  the  argument  of  perigee  combine  to  produce 
a  point  of  perigee  that  is  very  nearly  20  degrees  below  the 
descending  node  in  the  southern  hemisphere  of  the  earth's 
atmosphere.  The  results  of  the  sensitivity  analysis  plots 
for  the  right  ascension  of  the  ascending  node  and  the 
argument  of  perigee  are  shown  in  Figures  18  through  21. 

The  right  ascension  of  the  ascending  node  HO  is 
incremented  from  0  to  360  degrees.  Figure  18  presents  the 
10  kg  fuel  case  and  Figure  19  the  100  kg  fuel  case. 
Predicted  mission  life  ranges  are  ll  to  18  days  in  Figure 
18,  and  135  to  160  days  in  Figure  19. 

The  argument  of  perigee  GO  is  incremented  from  0  to  360 
degrees  with  a  10  kg  fuel  case  in  Figure  20,  and  a  100  kg 
fuel  case  in  Figure  21.  The  predicted  mission  life  are  13 
to  18  days,  and  148  to  160  days,  respectively. 

In  Figure  18  two  approximately  equal  mission  life  peaks 
and  valleys  are  observed.  The  valleys  correspond  to  the 
points  of  perigee  and  apogee  closely  approaching  the 
atmospheric  density  bulge.  The  peaks  occur  when  the  orbital 
plane  is  most  nearly  perpendicular  to  the  earth-sun  line 
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Mission  Life  versus  Right  Ascension 
of  the  Ascending  Node  (10  kg  of  Fuel) 
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Figure  19. 


Mission  Life  versus  Right  Ascension 
of  the  Ascending  Node  (100  kg  of  Fuel) 
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Figure  21. 


Mission  Life  versus  Argument 
of  Perigee  (100  kg  of  Fuel) 
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which  contains  the  atmospheric  density  bulge  center.  The 
change  in  shape  of  the  10  and  100  kilogram  curves  (Figures 
18  and  19)  is  due  to  the  movement  of  the  atmospheric  density 
bulge  during  mission  life  and  to  geopotential-induced 
orbital  drift. 

With  default  values  of  95  degrees  for  both  the  right 
ascension  of  the  ascending  node  and  inclination,  the  effect 
of  varying  the  argument  of  perigee  from  0  to  360  degrees  is 
to  move  the  point  of  perigee  around  the  initially  fixed 
orbital  plane  (see  Figure  4  in  Chapter  II)  .  The  plot  of 
mission  life  corresponding  to  this  motion  should  by  cyclic 
showing  one  peak  and  one  valley.  The  valley  corresponds  to 
the  perigee  point  closest  to  the  atmospheric  density-bulge 
center.  The  peak  in  plotted  mission  life  occurs  when  the 
apogee  point  is  cLosest  to  the  atmospheric  density -bulge 
center. 

Figures  20  and  21  show  the  plotted  results  of  the 
sensitivity  analysis  of  the  argument  of  perigee  GO.  The 
expected  single  peak  and  single  valley  are  seen.  The 
difference  in  the  10  kilogram  case  in  Figure  20  and  the  100 
kilogram  case  in  Figure  21  is  due  to  the  seasonal  movement 
of  the  atmospheric  density-bulge  and  geopotential-induced 
orbital  drift.  The  magnitude  of  the  mission  life  variations 
is  less  than  that  produced  by  varying  the  right  ascension  of 
the  ascending  node.  This  effect  is  due  to  low  eccentricity 
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orbit  (e  =0.01)  and  the  position  of  the  orbital  plane  with 
reference  to  the  earth-sun  line. 

Factors  which  primarily  move  the  orbital  plane  in 
relation  to  the  atmospheric  density-bulge  center  will  show 
two  peaks  and  two  valleys  in  their  single-cycle  mission-life 
range.  Factors  which  primarily  move  the  point  of  perigee  in 
an  initially  specified  orbital  plane  will  show  one  peak  and 
one  valley  in  their  single-cycle  range.  These 
considerations  assume  that  only  low  eccentricity  orbits  are 
being  considered. 

Based  on  the  above  discussion,  Figures  18  through  21 
appear  to  agree  with  expected  real-world  and  model- 
predicted  mission  life  trends. 

The  analysis  of  predicted  mission  life  trends  due  to 
varying  the  inclination  XI  is  more  complex.  Several  factors 
must  be  considered  in  addition  to  the  orbital-perigee 
atmospheric  density-bulge  encounter  geometry.  First,  the 
diurnal  atmosphere  model  used  by  the  propellant  longevity 
model  is  also  oblate,  and  is  thus  more  dense  at  the  equator 
than  at  the  poles.  This  fact  means  that  satellites  in  near 
polar  orbits  travel  through  an  overall  less  dense  atmosphere 
than  those  in  a  near  equatorial  orbit.  Second,  the 
atmosphere  rotates  in  the  same  direction  as  the  earth. 
Hence  prograde  orbital  satellites  experience  less  satellite- 
to-atmosphere  relative  velocity  than  satellites  in 
retrograde  orbits.  If  the  satellite  altitude  is  low  enough, 


a  near-equatorial  prograde  orbit  satellite  has  a  longer 
mission  life  than  the  same  satellite  in  a  near-equatorial 
retrograde  orbit.  These  additional  effects  complicate  the 
analysis  of  mission  life  resulting  from  varying  the  orbital 
inclination. 

The  default  value  position  of  the  center  of  the 
atmospheric  density-bulge  is  300  degrees  longitude  and  20 
degrees  south  latitude.  The  bulge  is  moving  at  a  rate  of 
about  one  degree  longitude  per  day  toward  the  vernal 
equinox.  The  point  of  perigee  is  default  positioned  at 
approximately  275  degrees  longitude  and  20  degrees  south 
latitude.  As  can  be  seen,  the  point  of  perigee  is  quite 
close  to  the  center  of  the  atmospheric  density-bulge.  The 
motion  of  the  point  of  perigee,  as  the  inclination  is 
increased  from  near  0  degrees  to  near  180  degrees,  describes 
a  semicircle  20  degrees  in  radius  around  the  descending 
node.  The  point  of  perigee  remains  quite  close  to  the 
initial  position  of  the  center  of  the  atmospheric  density- 
bulge. 

Figures  22,  23,  24,  and  25  are  the  predictions  of  the 
sensitivity  analysis  plots  of  orbital  inclination  XI.  The 
incremented  range  is  from  22.5  to  157.5  degrees.  Nearer 
equatorial  inclinations  are  prohibited  by  the  available 
arithmetic  precision  on  the  IBM  3033.  Figure  22  shows  the 
10  kg  fuel  case,  Figure  23  the  100  kg  fuel  case,  and  Figure 
24  the  100  kg  fuel  case  with  a  granularity  of  1/100  and  a 
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rotation  factor  of  0.9.  Predicted  mission  life  ranges  are 
7.8  to  13.7  days  in  Figure  22,  and  90  to  138  days  in  Figures 
23  and  25. 

Figure  25  shows  the  effect  of  the  using  100  kilograms  of 
compensating  fuel  and  a  rotation  factor  of  0.9.  This 
situation  is  compared  to  100  kilograms  of  compensating  fuel 
and  a  rotation  factor  of  1.0  in  Figure  23.  Only  very  slight 
differences  occur  between  Figures  23  and  25.  Mission  life 
is  a  little  longer  in  the  retrograde  orbits  using  a  rotation 
factor  of  0.9  and  a  little  shorter  in  the  prograde  orbits 
using  a  rotation  factor  of  0.9. 

When  the  satellite  orbit  is  lowered  to  an  altitude  at 
perigee  of  200  kilometers  (AO  =  1.03)  with  1000  kg  of  fuel, 
mission  life  is  reduced  in  the  retrograde-equatorial  orbit 
below  that  of  the  prograde-equatorial  orbit.  The  observed 
peaks  at  70  and  140  degrees  also  shift  leftward.  The 
observed  dimple  seen  in  Figure  22  at  the  critical 
inclination  of  63  degrees  may  be  geopotential-induced  or  an 
artifact  of  the  Kozai  approximations  as  discussed  by  Taff 
[Ref.  2:pp  332-342]. 

A  further  explanation  of  the  plot  of  mission  life  versus 
inclination  shown  in  Figure  22  is  presented  here.  A  more 
detailed  study  of  this  input  variable  is  worthy  of  a 
complete  paper  itself. 

The  low  mission  life  found  at  low  values  of  inclination 
is  due  to  the  combination  of  the  point  of  perigee  position 


and  oblate-atmosphere  effects.  The  low  inclination  orbits 
travel  through  the  "fat  part"  of  the  oblate  atmosphere  near 
the  equator.  As  the  inclination  increases,  the  point  of 
perigee  moves  away  from  the  center  of  the  atmospheric 
density-bulge  and  the  satellite's  orbit  moves  up  out  of  the 
thick  equatorial  atmosphere.  These  two  factors  combine  to 
produce  the  observed  rapid  rise  in  mission  life  seen  in 
Figure  22.  Beyond  near-polar  values  of  the  inclination,  the 
predicted  mission  life  again  falls.  The  mission  life  valley 
seen  at  120  degrees  is  probably  due  to  maximum  orbital  time 
near  the  atmospheric  density-bulge  center. 

Figure  23  shows  the  effect  of  starting  with  100 
kilograms  of  initial  drag-compensating  fuel.  Note  the  many 
peaks  and  valleys.  To  confirm  that  these  peaks  and  valleys 
represent  actual  predicted  mission  life  trends  and  not 
simply  data  scatter,  the  granularity  was  decreased  to  1/100. 
Figure  24  is  a  plot  of  these  results,  and  confirms  the 
trends  outlined  in  Figure  23.  The  many  peaks  and  valleys 
superimposed  over  the  general  trend  of  mission  life  seen  in 
Figure  22  are  probably  due  to  geopotential-induced  orbital 
drift. 

The  results  of  the  sensitivity  and  trend  analysis  of  the 
orbital  inclination  XI  have  been  discussed  with  reference  to 
in  Figures  22  through  25.  Future  investigation  of  the 
plotted  data  may  reveal  more  insight  toward  understanding. 


the  underlying  mechanisms.  The  sensitivity  analysis  plots 
show  no  trends  that  appear  unreasonable. 

The  sensitivity  analysis  for  the  input  variables 
effecting  the  orbital-perigee  atmospheric  density-bulge 
encounter  geometry  (category  one)  have  now  been  presented. 
The  next  set  of  analyses  are  those  of  the  second  category, 
which  contains  the  input  variables  affecting  the  slowing 
power  of  the  atmosphere.  These  variables  are:  decimetric 
solar  flux  F107,  drag  coefficient  CD,  cross-sectional  area 
A,  specific  impulse  SPI,  and  the  atmospheric  rotation  factor 
D.  The  initial  mass  of  drag-compensating  fuel  is  also 
included  in  this  category  because  it  improves  satellite 
mission-life  resistance  to  atmospheric  drag. 

The  real-world  expected  effect  on  mission  life  as  a 
result  of  varying  the  decimetric  solar  flux  F107  is  to  vary 
the  heating  effectiveness  of  the  sun.  As  the  decimetric 
solar  flux  increases,  the  density  of  the  atmosphere  will 
increase . 

The  Jacchia  J60  model  atmosphere  used  in  the  propellant 
longevity  model  subroutine  JAC60  uses  the  F107  flux  as  a 
multiplier  for  the  exponential  density  term.  The  expected 
model  predictions  of  mission  life  as  the  F107  flux  is  varied 
through  its  range  of  50  to  300  is  a  multiplicatively 
decreasing  plot. 

The  results  of  the  sensitivity  analysis  for  predicted 
mission  life  as  the  decimetric  solar  flux  is  varied  are 


iiwuiBimiin 


presented  in  Figures  26  and  27.  Figure  26  is  the  result  of 
10  kg  of  initial  maneuvering  fuel,  and  Figure  27  is  the 
result  with  100  kg  of  fuel.  The  predicted  mission  life  in 
Figure  26  ranges  from  24.99  to  3.96  days.  The  100  kg  fuel 
case  in  Figure  27  predicts  a  mission  life  ranging  from 
270.12  to  43.07  days.  These  results  are  in  agreement  with 
expectations . 

The  next  set  of  analyses  consider  effects  of  the  drag 
coefficient,  satellite  cross-sectional  area,  and  specific 
impulse.  All  of  these  factors  are  multipliers  acting  to 
change  the  resistance  of  the  satellite  to  atmospheric  drag. 
The  drag  coefficient  CD  is  a  proportionality  constant  in  the 
drag  force  model.  The  satellite  cross-sectional  area  A 
changes  the  size  of  the  satellite  that  the  atmosphere 
actually  acts  against.  As  satellite  cross-sectional  area 
increases  it  will  take  multiplicatively  more  energy  to 
compensate  for  the  orbital  energy  lost  by  the  satellite 
working  against  the  drag  force.  The  specific  impulse  SPI  is 
a  multiplier  that  expresses  the  effectiveness  of  the 
expelled  fuel  mass  in  overcoming  the  drag  force. 

The  rotation  factor  D  is  the  ratio  of  the  atmospheric 
rotation  rate  to  the  earth's  rotation  rate.  Its  effect  is 
to  change  the  satellite  to  atmospheric  relative  velocity  as 
the  inclination  changes.  A  logical  range  of  D  is  between 
the  values  0.9  and  1.0.  The  sensitivity  analysis  range  of 
0.1  to  10.0  was  chosen  to  show  model  sensitivity  for  a  wide 
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range  of  D.  As  presented  in  Chapter  II  Equation  12,  the  D 
factor  is  combined  with  the  inclination  and  the  ratio  of  the 
radius-at-perigee  to  the  velocity-at-perigee  producing  the 
6  factor.  This  factor  is  used  to  convert  the  satellite 
tangential  velocity  into  the  satellite- to- atmospheric 
relative  velocity. 

The  propellant  longevity  model  uses  A,  CD,  <$  ,  and  SPI 
as  proportionality  constants,  multiplicatively  changing  the 
effectiveness  of  the  expelled  fuel  mass.  The  routine  PLEP 
combines  these  four  factors  into  a  single  constant  DEL  given 
by: 

DEL  -  A  *  CD  *  6  *  3 . 14159/SPI  (14) 

This  constant  DEL  is  passed  to  the  routine  KOZAI  which  in 
turn  passes  it  to  the  routine  DRAG.  In  the  routine  DRAG, 
DEL  is  used  as  the  final  multiplier  in  determining  required 
fuel  for  each  orbit. 

The  expected  model-predicted  mission-life  as  a  result  of 
changing  these  factors  is  a  multiplicative  increase  when  A, 
CD,  or  <5  are  increased  linearly.  A  multiplicative  decrease 
should  occur  when  SPI  is  increased  linearly. 

Expected  results  of  varying  D  from  0.1  to  10.0  in  the 
near-polar  retrograde  orbit  are  shown  in  Figure  28.  Note 
that  a  very  slight  lowering  of  mission  life  is  seen  when  the 
rotation  factor  increases.  Figures  29  and  30  show  the 


plotted  result  of  varying  D  through  the  same  range  of  0.1  to 
10.0.  However,  in  Figure  29  an  inclination  of  25  degrees  is 
used  instead  of  the  default  95  degrees,  and  in  Figure  30  an 
inclination  of  155  degrees  is  used.  In  Figure  29  it  is  seen 
that  the  mission  life  multiplicatively  increases  as  D  is 
linearly  increased.  Figure  30  shows  a  multiplicative 
decrease  in  mission  life  as  D  is  linearly  increased.  These 
trends  are  due  to  the  change  in  the  satellite-to-atmosphere 
relative  velocity.  The  increasing  mission  lif~  expectancy 
within  a  strongly  prograde  orbit  (Figure  29)  is  a  result  of 
the  decreasing  satellite-to-atmosphere  relative  velocity. 
The  decreasing  mission  life  expectancy  within  a  strongly 
retrograde  orbit  (Figure  30)  is  a  result  of  increasing 
satellite-to-atmosphere  relative  velocity.  The  plots  shown 
in  Figures  28,  29  and  30  are  in  agreement  with  expected 
results . 

The  results  of  the  sensitivity  analysis  plots  on  mission 
life  as  a  result  of  varying  CD,  A,  and  SPI  are  shown  in 
Figures  31  through  35. 

Figures  31  and  32  result  from  varying  the  satellite 
drag  coefficient  CD  from  1.0  to  3.5.  Figure  31  represents 
the  10  kg  fuel  case  and  Figure  32  the  100  kg  case.  The 
predicted  mission  life  range  in  Figure  31  is  31.18  to  8.63 
days,  and  in  Figure  32  from  308  to  99.58  days. 


The  sensitivity  analysis  mission-life  trend  plots  of 
satellite  cross-sectional  area  A  are  shown  in  Figures  33  and 
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Figure  34.  Mission  Life  versus  Cross-Sectional 
Area  (100  kg  of  Fuel) 
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Figure  35.  Mission  Life  versus  Specific 
Impulse  (10  kg  of  Fuel) 


34.  Figure  33  is  the  10  kg  fuel  case  and  Figure  34  the  100 
kg  case.  In  both  cases,  A  was  varied  from  0.000010  to 
0.0001000  square  kilometers.  The  ranges  mapped  by  these 
plots  are  80  to  0.64  days  and  736.29  to  108.7  days, 
respectively . 

Figure  35  is  the  analysis  of  motor  specific  impulse  SPI 
varying  from  200  to  3000  with  10  kg  of  initial  fuel.  The 
apparent  linear  shape  of  the  SPI  plot  is  due  to  its  inverse 
multiplicative  action  and  the  range  of  the  constant  DEL  (14) 
over  which  it  acts.  As  SPI  increases,  the  drag  compensating 
fuel  is  more  effective.  Linearly  increasing  SPI  has  the 
same  effect  as  starting  at  the  default  value  *  on  the  cross- 
sectional  plots  and  moving  along  the  plot  to  the  left.  Note 
that  the  shape  of  cross-sectional  mission  life  plot  in  this 
range  is  closer  to  the  expected  mirror  image  of  the  SPI 
plot. 

As  noted  earlier,  mission  life  plots  showing  "knees” 
have  significance  to  research  and  design  studies.  The  knees 
in  the  cross-sectional  area  mission  life  plots  are  very 
pronounced.  Less  distinctive  knees  are  seen  in  the  mission- 
life  plots  of  the  decimetric  solar  flux  and  the  drag 
coefficient.  The  default  value  of  the  cross-sectional  area 
is  right  in  the  knee  range  of  the  predicted  mission-life 
plots,  as  seen  in  Figures  34  and  35. 

The  final  input  variable  in  the  second  category  is  the 
initial  propellant  mass  W.  The  expected  mission  life  of  a 
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drag  compensating  satellite  as  fuel  is  increased  linearly 
should  be  a  linear  increase  in  the  predicted  mission  life. 
A  given  amount  of  fuel  produces  a  certain  mission  life. 
Twice  that  amount  of  fuel  in  the  same  general  configuration 
should  produce  about  twice  the  mission  life.  The  expected 
shape  of  the  sensitivity  analysis  plot  of  predicted  mission 
life  for  linearly  increasing  initial  propellant  mass  is  a 
approximately  straight  line.  Figure  36  is  the  result  of 
varying  W  on  mission  life.  The  plot  is  approximates  a 
straight  line  as  expected.  The  ripples  are  probably  due  to 
geopotential-induced  orbital-drift  or  to  movement  of  the 
atmospheric  density-bulge  during  mission  life. 

The  results  of  the  sensitivity  and  trend  analysis 
mission-life  plots  for  CD,  A,  D,  SPI,  and  W  agree  with 
expectations.  All  input  variables  in  the  second  category 
have  been  examined.  The  next  input  variables  to  be  examined 
are  in  the  third  category. 

Consider  now  the  input  variables  that  affecting  the 
depth  of  penetration  of  the  satellite  into  the  atmosphere  at 
perigee.  The  variables  that  affect  the  satellite  altitude 
at  perigee  are  the  semi-major  axis  and  the  orbital 
eccentricity.  These  variables  show  the  most  pronounced 
effects  on  mission  life.  The  real-world  effect  on  mission 
life  when  linearly  changing  the  altitude  at  perigee  is  an 
exponential  increase  in  mission  life.  The  effect  of  varying 
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the  eccentricity  and  the  semi-major  axis  on  the  altitude  at 
perigee  is  summarized  by 

rp  =  a  *  (1  -  e)  (15) 

Here,  rp  is  the  radius  at  perigee,  a  the  semi-major  axis, 
and  e  the  eccentricity.  A  linear  increase  in  the  semi¬ 
major  axis  results  in  a  modified  exponential  increase  in 
mission  life.  A  linear  increase  of  the  eccentricity  results 
in  a  modified  exponential  decrease  in  mission  life. 

The  propellant  longevity  model  computes  the  satellite 
altitude  at  perigee  in  the  subroutine  SALT.  The  subroutine 
DRAG  then  passes  the  altitude  at  perigee  to  the  routine 
JAC60  where  it  is  used  to  develop  an  exponent  for  an 
empirical  model  of  atmospheric  density. 

Mission  life  results  of  varying  the  orbital  eccentricity 
ES  from  0.01  to  0.1  are  presented  in  Figures  37  and  38. 
Figure  37  is  the  10  kg  fuel  case,  and  Figure  38  the  100  kg 
fuel  case.  The  predicted  mission  life  ranges  from  14  to  0 
days  in  Figure  37,  and  150  to  0  days  in  Figure  38. 

An  analysis  of  the  semi-major  axis  A0  is  presented  in 
Figure  39  with  a  range  of  1.025  to  1.175  earth  radii.  In 
this  case,  1  kg  of  fuel  is  allocated  initially.  The 
resulting  mission  life  is  from  0  to  384  days. 

The  plots  of  mission  life  resulting  from  varying  the 
semi-major  axis  A0  and  the  orbital  eccentricity  ES,  show 
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Figure  37.  Mission  Life  versus  Orbit 

Eccentricity  (10  kg  of  Fuel) 


major  knees  in  their  curves.  The  *  symbol  indicates  the 
default  value  of  the  varied  element.  Note  that  the  amount 
of  fuel  used  in  the  analysis  of  the  semi-major  axis  is  l  kg. 
The  use  of  the  usual  10  kg  produces  a  plot  that  exceeding 
the  reasonable  range  of  mission  life.  The  results  shown  in 
Figures  37  through  39  are  in  agreement  with  expectations. 

The  input  variables  in  the  fourth  category  (no  effect  on 
mission  life)  are  the  mean  anomaly  ETU,  initial  satellite 
weight  WT,  average  decimetric  solar  flux  FBAR,  and  the 
geomagnetic  index  AKP.  FBAR  and  AKP  have  no  effect  on 
mission  life  due  to  the  limitations  of  the  Jacchia  J60 
atmosphere  model.  In  the  real  world  these  two  factors  may 
combine  with  the  decimetric  solar  flux  F107  to  produce 
changes  in  the  atmospheric  density  of  three  orders  of 
magnitude. 

The  model  theory  used  by  Dr.  Parks  [Ref.  1]  to  predict 
the  propellant  longevity  of  a  satellite  doing  intrack  micro- 
thrusting  to  overcome  atmospheric  drag  is  independent  of  the 
satellite  mass.  The  fuel  required  to  overcome  drag-induced 
orbital  energy  loss  is  based  strictly  on  the  satellite's 
shape  and  size. 

Figure  40  presents  the  sensitivity  analysis  plot  of 
predicted  mission  life  as  a  result  of  varying  the  satellite 
initial  mass  WT.  The  plot  is  a  straight  horizontal  line,  as 
expected . 
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Figure  40.  Mission  Life  versus  Satellite 
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The  final  analysis  is  concerned  with  the  mean  anomaly. 
In  the  near-circular  default  orbit,  the  mean  anomaly  and  the 
eccentric  anomaly  are  nearly  the  same.  The  Keplerian  mean 
anomaly  functions  to  position  the  satellite  in  its 
designated  orbit.  For  a  mission  life  lasting  several  days, 
it  should  make  no  difference  where  the  satellite  is 
initially  positioned  in  its  orbit.  Rather,  it's  the  orbit 
itself  that  makes  the  difference. 

Analysis  plots  of  the  m.ean  anomaly  ETU  are  presented  in 
Figures  41  and  42.  The  incremented  range  is  from  0  to  179 
degrees.  Figure  41  is  the  10  kg  fuel  case  and  Figure  42  is 
the  100  kg  fuel  case.  Predicted  results  range  from  13.74  to 
13.87  days  in  Figure  41,  and  from  149  to  154  days  in  Figure 
42. 

The  results  in  Figure  41  are  as  expected,  an  approximate 
horizontal  line.  In  Figure  42,  a  valley  is  seen  in  the  40 
to  120  degree  range.  The  valley  is  only  four  days  deep,  but 
is  difficult  to  explain.  Further  research  may  reveal  its 
underlying  mechanism. 

The  results  of  the  trend  and  sensitivity  analysis  have 
been  presented.  Mission  life  has  been  found  to  be  most 
sensitive  to  those  factors  that  affecting  the  orbital  radius 
at  perigee.  The  main  unexpected  results  are  found  in  the 
analysis  of  the  inclination  and  the  mean  anomaly.  The 
pattern  of  many  peaks  and  valleys  found  in  the  100  kilogram 
fuel  case  (Figures  23,  24  and  25)  warrants  further  research. 


MEAN  ANOMALY  (DEGREES) 


Figure  42. 


Mission  Life  versus  Mean 
Anomaly  (100  kg  of  Fuel) 
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The  four  day  dimple  seen  in  the  mean  anomaly  plot  (Figure 
41)  also  warrants  further  study. 


V.  SUMMARY 


This  final  chapter  presents  the  summary  and  conclusions 
of  the  thesis.  Recommendations  for  model  use  are  also 
given . 

The  low-earth-orbit  satellite  propellant  longevity  model 
has  been  tested  in  a  sensitivity  analysis.  Based  on  the 
results  of  the  sensitivity  analysis  the  model  reasonableness 
is  confirmed.  The  first  chapter  presented  a  model  overview, 
model  motivation,  and  the  process  employed  to  verify  the 
model.  Chapter  I  also  discusses  some  changes  made  to  the 
model  computer  program  for  the  IBM  3033  at  the  Naval 
Postgraduate  School.  Chapter  II  presents  the  orbital 
mechanics  and  atmospheric  concepts  as  background  for  model 
understanding  and  it's  operation.  In  Chapter  III  the  model 
and  its  computer  program  were  discussed.  Operating 
instructions  specific  for  the  computer  program  use  at  the 
Naval  Postgraduate  School  are  also  presented  there.  The 
sensitivity  analysis  results  are  presented  in  Chapter  IV. 

The  propellant  longevity  model  predicts  how  long  the 
fuel  of  a  low-earth-orbital  satellite  thrusting  to  overcome 
atmospheric  drag  will  last  before,  starting  uncontrolled 
decay  into  the  atmosphere.  Model  verification  reasonable¬ 
ness  tests  reveal  no  inconsistencies  with  the  established 
theory.  The  results  of  the  sensitivity  analysis  are 
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therefore  in  agreement  with  the  expectations.  The 
irregularities  in  the  model  predictions  (occurring  with 
inclinations  of  less  than  22.5  degrees)  are  due  to  computer 
round  off  error,  not  the  model  itself. 

The  propellant  longevity  model  currently  uses  the 
Jacchia  J60  model  atmosphere.  This  atmosphere  model  may 
significantly  reduce  the  real-world  accuracy  of  predicted 
mission  lifetimes. 

The  short  term  confidence  factor  in  predicted  mission 
life  might  be  increased  by  the  employment  of  the  Jacchia  J77 
model  atmosphere.  One  drawback  to  the  use  of  the  Jacchia 
J77  model  atmosphere  is  that  it  may  double  the  computer  run 
time. 

One  solution  to  the  trade-off  of  increased  accuracy  of 
the  Jacchia  J77  model  versus  the  shorter  computation  time  of 
the  Jacchia  J60  model,  would  be  to  allow  the  user  to  select 
an  atmosphere  option. 

Transportation  of  the  computer  program  of  the  model  to 
an  IBM  AT  could  provide  more  flexibility  in  preliminary 
design  analysis. 

This  model  is  a  tool  for  the  low-earth-orbit  satellite 
system  planner  in  designing  and  managing  satellite  systems 
that  have  precise  orbital  element  control  requirements. 
These  systems  are  those  with  SDI  and  reconnaissance 


This  model,  with  the  addition  of  the  Jacchia  J77  model 
atmosphere  and  a  subroutine  to  calculate  fuel  required  to 
compensate  for  aspherical  geopotential  forces,  could  be  used 
to  optimize  SDI  defensive  satellite  constellation  design. 

Existing  real-world  data  to  compare  with  the  model 
predictions  is  limited  to  a  very  classified  system.  The 
final  step  in  verifying  the  propellant  longevity  model  would 
be  to  make  this  comparison  test. 
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APPENDIX  A 

COMPUTER  PROGRAM  LISTING 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


XX  XX 
XX  PROPELLENT  LONGEVITY  xx 
XX  FOR  xx 
XX  LOW-ALTITUDE  EARTH  ORBIT  SATELLITES  xx 
xx  (A.D.  PARKS)  xx 

.  XX  XX 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  X 

X  THIS  PROGRAM  IS  THE  CODING  OF  THE  THEORY  CONTAINED  IN  NSWC  X 

X  TR  83-243.  IT  WAS  USED  TO  GENERATE  THE  NUMERICAL  EXAMPLES  ON  X 

X  PAGES  13  TO  17.  INITIALLY  CODED  FOR  THE  CDC  6700  IN  FORTRAN  x 

X  IV  BY  A.D.  PARKS,  THE  CODE  HAS  BEEN  TRANSPORTED  TO  THE  IBM  3033  x 

X  AND  TRANSLATED  TO  FORTRAN  77  BY  LT.  C.D.  NOBLE.  ADDITIONAL  DOC-  x 

X  UMENTATION  WAS  ADDED  DURING  THE  TRANSPORTATION  AND  TRANSLATION.  x 

x  x 

X  THE  REFERENCE  MATERIAL  CITED  IN  THE  CODE  IS  MOSTLY  FORM  THE  x 

X  NAVAL  SURFACE  WEAPONS  CENTER  (NSWC)  TECHNICAL  REPORTS  (TR) .  x 

X  NSWC  TR  *  ARE  CITED  IN  THE  CODE  AS  COMMENTS.  OTHER  REFERENCES  x 

X  ARE  SITED  BY  NUMBER.  THEY  ARE.  x 

X 

BROUWER,  D.,  "SOLUTION  OF  THE  PROBLEM  OF  x 
ARTIFICIAL  SATELLITE  THEORY  WITHOUT  DRAG",  X 
THE  ASTRONOMICAL  JOURNAL,  VOL  64,  NO  1274,  x 
PP.  378-397,  1959.  « 

x 

AGRAWAL,  B.N.,  "DESIGN  OF  GEOSYNCHRONOUS  x 

SPACECRAFT",  PRENTICE-HALL,  1986.  x 

x 

JACCHIA,  L.G.,  "THERMOSPHERIC  TEMPERATURE,  X 
DENSITY  AND  COMPOSITION.  NEW  MODELS",  x 

SPECIAL  REPORT  375,  SMITHSONIAN  ASTRO-  x 

PHYSICAL  OBSERVATORY,  1977.  x 

x 

JACCHIA,  L.G.,  "EMPIRICAL  MODELS  OF  THE  x 

THERMOSPHERE  AND  REQUIREMENTS  FOR  IMP-  x 

ROVEMENTS",  AFGL-TR-81-0195,  ADA101946,  x 

1981 . 

x 

JACCHIA,  L.G.,  "ANALYSIS  OF  DATA  FOR  THE  x 
DEVELOPMENT  OF  DENSITY  AND  COMPOSITION  x 

MODELS  OF  THE  UPPER  ATMOSPHERE",  AFGL-TR-  x 
81-0230,  ADI 06420 ,  1981.  x 

x 

ORR,  L.H.,  "ORBITAL  LIFETIME  PROGRAM",  x 

NASA  TM  87587,  1985.  x 

x 

TAFF,  L.G.,  "CELESTIAL  MECHANICS.  A  COMP-  x 
UTATIONAL  GUIDE  FOR  THE  PRACT I  ONER" ,  x 

JOHN  WILEY  S  SONS.  INC.,  1985.  x 

x 

FLEAGLE,  R.G.,  "AN  INTRODUCTION  TO  ATMOS-  x 

PHERIC  PHYSICS",  ACADEMIC  PRESS,  1980.  x 
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REFERENCE  1 

REFERENCE  2 
REFERENCE  3 
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*  MODULE 
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X 


X 
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X 
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X 
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X 

X 

X 

X 
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X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


PLEP 


POSVEL 


NWTRPH 


FUNCTIONS 

XSPA-XSPR 


MEAN 


GEOP 


SOLOC 


FUNCTION 


SYSTEM  INITIALIZATION  AND  CONTROL, 


CALLED  BY. 


CALLS. 


COMPUTES  SATELLITES  POSITION  AND  VELOCITY  IN 
INERTIAL  SPACE  FROM  THE  OSCULATING  ORBITAL 
ELEMENTS.  NSHC  TR  82-387  PP.  24-25. 


CALLED  BY. 


SATLOC 


CALLS. 
(SUBROUTINES) 
NWTRPH 


TIMES  CALLED 
1 


A  NEWTON-RAPHSON  SOUUTION  TO  KEPLERS  EQUATION. 
NSWC  TR  82  387  P  23. 

CALLS. 

(SUBROUTINES)  TIMES  CALLED 
xxxx  xxxx 


CALLED  BY. 


POSVEL 


COMPUTES  SHORT  PERIODIC  EFFECTS  ON  ORBITAL 
ELEMENTS.  CALLED  BY  SUBROUTINE  MEAN.  NSWC  81-456 
PAGES  9-11  AND  NSWC  TR  83  243  REF  3. 

CALLS. 

(SUBROUTINES) 
xxxx 


CALLED  BY. 


MEAN 


TIMES  CALLED 
xxxx 


USES  THE  WALTER  METHOD  TO  CONVERT  OSCULATING 
ORBITAL  ELEMENTS  TO  KOZAI  LIKE  MEAN  ELEMENTS 
USING  FIRST  ORDER  PERIODIC  VARIATIONS. 

SIMILAR  TO  THE  METHOD  IN  NSWC  TR  83  135  PP  16-18. 
CALLED  BY:  CALLS: 

(FUNCTIONS) 

KOZAI  XSPA-XSPR 


TIMES  CALLED 
1 


COMPUTES  AVERAGED  SECULAR  CHANGES  IN  KOZAI  MEAN 
ELEMENTS  INDUCED  BY  GRAVITATIONAL  FIELD  THROUGH 
FOURTH  ZONAL  HARMONIC  CJ2-J4  TERMS).  NSWC  TR  81- 
456.  EQS.C2.5). 

CALLED  BY.  CALLS. 

(SUBROUTINES)  TIMES  CALLED 
KOZAI  xxxx  xxxx 


COMPUTES  THE  RIGHT  ASCENSION  AND  THE  DECLINATION 
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X 

(SUBROUTINES) 

TIMES  CALLED 

X 

X 

xxxx 

SNALYS 

1 

X 

X 

x 

KOZAI 

1/ITERATION 

X 

X 

SNALYS 

SENSITIVITY  MODE  SCANNED  ELEMENT 

GENERATOR. 

X 

X 

X 

CALLED  BY. 

CALLS. 

X 

X 

(SUBROUTINES) 

TIMES  CALLED 

X 

X 

£ 

PLEP 

xxxx 

xxxx 

X 

X 

KOZAI 

SYSTEM  WORKHORSE, 

DRIVER  AND  RECORD  KEEPER. 

X 

X 

X 

NSWC  TR  83-135. 

X 

X 

CALLED  BY. 

CALLS. 

X 

X 

(SUBROUTINES) 

TIMES  CALLED 

X 

X 

PLEP 

BRAUER 

1 

X 

X 

PERIOD 

1 

X 

X 

MEAN 

1 

X 

X 

THRST 

1/ ( SAT . REV . ) 

X 

X 

. 

GEOP 

1/ (SAT. REV. ) 

X 

X 

jf 

DRAG 

1/($AT . REV . ) 

X 

X 

BRAUER 

CONVERTS  NAVSPASUR 

BROUWER  MEAN  ORBITAL  ELEMENT 

X 

X 

X 

SET  INTO  A  BROUWER 

OSCULATING  ORBITAL  SET. 

X 

X 

A  CODING  OF  THE  BROUWER-LYDDANE  THEORY. 

X 

X 

£ 

NSWC  TR  83-107,  83 

-135,  BROUWER  (REF.  1). 

X 

X 

X 

CALLED  BY. 

CALLS. 

w 

X 

X 

(SUBROUTINES) 

TIMES  CALLED 

X 

X 

KOZAI 

xxxx 

xxxx 

X 

x 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 
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X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


PERIOD 


DRAG 


JAC60 


THRST 


SALT 


SATLOC 


FUNCTION 

(MMBSIR) 


OF  THE  SUN  FOR  THE  MODEL  ATMOSPHERE. 
456  PP.  20-21. 

CALLED  BY i  CALLS: 

(SUBROUTINES) 
DRAG  xxxx 


NSWC  TR  81- 


TIMES  CALLED 
xxxx 


COMPUTES  THE  ANOMALISTIC  PERIOD. 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
KOZAI  xxxx 


TIMES  CALLED 
xxxx 


CALLED 


COMPUTES  DRAG  EFFECTS  AND  FUEL  MASS  DECREMENT 
REQUIRED  TO  MAINTAIN  THE  SEMIMAJOR  AXIS  VIA 
IN  TRACK  MICROTHRUSTING .  NSWC  TR  82-243  EQ . 

58.  THIS  ONE  EQUATION  IS  THE  HEART  OF  THE  MODEL. 
CALLED  BY:  CALLS: 

(SUBROUTINES) 

KOZAI  SOLOC 

SATLOC 
SALT 
JAC60 

(FUNCTIONS) 

MMBSIR 
XSPM 


TIMES 

1 

3 

1 

3 


COMPUTES  ATMOSPHERIC  DENSITY  BASED  ON  THE  JACCHIA 
(60)  MODEL  ATMOSPHERE.  CALLED  THREE  TIMES  TO 
COMPUTE  DIFFERENT  VALUES  FOR  THE  DENSITY  FOR  EACH 
REVOLUTION:  MAX,  MIN  AND  OP. 


CALLED  BY: 


DRAG 


CALLS: 
(SUBROUTINES) 
xxxx 


TIMES  CALLED 
xxxx 


CALCULATES 
AN  ORBITAL 


CHANGES  TO  KOZAI  ORBITAL  ELEMENTS  IF 
ADJUST  IS  SCHEDULED.  CALCULATES  THE 
MASS  OF  FUEL  USED  TO  PERFORM  THE  ADJUST. 

NSWC  TR  83-31  EQ.  16,  AND  NSWC  TR  81-456 
CALLED  BY.  CALLS-. 

(SUBROUTINES) 

KOZAI  xxxx 


EQ.  2-26 


TIMES  CALLED 
xxxx 


COMPUTES  SATELLITE  ALTITUDE.  NSWC  TR  81-456  EQ 
4.34. 

CALLED  BY.  CALLS: 

(SUBROUTINES) 

DRAG  xxxx 


TIMES  CALLED 
xxxx 


COMPUTES  SATELLITE  LOCATION.  RIGHT  ASCENSION, 
DECLINATION  AND  ITS  GEOMAGNETIC  LATITUDE. 

CALLED  BY:  CALLS: 

(SUBROUTINES)  TIMES  CALLED 
DRAG  POSVEL  1 


COMPUTES  BESSEL  FUNCTIONS.  AN  IM3L  ROUTINE 
RESIDENT  ON  THE  IBM  3800.  THIS  HILL  VARY  WITH 
LOCAL  INSTALLATIONS. 

CALLED  BY:  CALLS: 

(SUBROUTINES)  TIMES  CALLED 
DRAG  xxxx  xxxx 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

**  PROGRAM  PLEP  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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PROGRAM  PLEP 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X  THIS  MODULE  IS  THE  INPUT,  CONTROL  AND  INITIALIZATION  SECTION  FOR  THE* 

*  REST  OF  THE  SYSTEM.  IT  ACCEPTS  THE  'FIVE  CARD'  FORMAT  OF  NAVSPASUR  * 

*  AS  THE  BROUWER  MEAN  ORIBITAL  ELEMENT  SET.  THE  SPECIFIC  SATELLITE  * 
X  PARAMETERS  ARE  ENTERED.  ATMOSPHERIC  MODEL  ELEMENTS  ARE  INITIALIZED.  * 

*  FINALLY,  SYSTEM  CONTROL  CONSTANTS  ARE  SET  TO  GENERATE  THE  OUTPUT.  A  X 

X  FURTHER  DESCRIPTION  OF  THE  ABOVE  FUNCTIONS  IS  CONTAINED  IN  THE  * 

X  DESCRIPTION  OF  ARGUMENTS.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  X 

*  DESCRIPTION  OF  ARGUMENTS  x 


INPUT  ARGUMENTS 
PROGRAM  CONTROLS 


ITERl. 

ITER2 


PSENS 


IRV(K). 


DA(K) 


PSIZE. 


INTEGER,  MAXIMUM  NUMBER  OF  REVOLUTIONS  DESIRED  * 
FOR  CONSIDERATION.  A  CONTROL  PARAMETER.  * 

INTEGER,  ORBITAL  PRINT  FREQUENCY .  A  CONTROL  * 

"PARAMATER.  SETTING  ITER2  =  TO  100  WILL  REDUCE  * 
TABULAR  OUTPUT  BY  ONLY  PRINTING  OUTPUT  EVERY  100  x 
ORBITS.  x 

INTEGER,  CONTROL  PARAMETER.  PSENS=1  SPECIFIES  X 
'SENSITIVITY  ANALYAIS  MODE.  PSENS=0  SPECIFIES  * 

STANDARD  PROPELLENT  LONGIVITY  MODE.  * 

INTEGER.  CONTROL  PARAMETER.  IET=1  PRINTS  THE  * 

'INITIAL  ANOMALISTIC  PERIOD.  IET*0  SUPRESSES  PRINT.* 
SYSTEM  RESET  TO  0  WHEN  SENSITIVITY  MODE  IS  SELECTED 
.INTEGER,  NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS.  THIS  x 
IS  A  CONTROL  PARAMETER  WITH  ALLOWED  VALUES  FROM  * 
ZERO  TO  TEN.  THE  ONLY  ORBIT  ADJUSTS  CURRENTLY  * 
SUPPORTED  ARE  CHANGES  TO  THE  SEMI-MAJOR  AXIS.  * 

.INTEGER  ARRAY.  (K  =  1Q  MAX). THE  ELEMENTS  OF  THE  x 
ARRAY  SPECIFY  THE  REVOLUTION  NUMBER  DURING  WHICH  * 

EACH  OF  THE  ORBITAL  ADJUSTS  ARE  PERFORMED.  K  MUST  * 
EQUAL  NOA.  * 

INTEGER  ARRAY.  CK*10  MAX)  THE  ELEMENTS  OF  THE  x 

ARRAY  SPECIFY  THE  CHANGE  IN  KILOMETERS  OF  THE  * 

SEMIMAJOR  AXIS.  K  MUST  EQUAL  NOA.  * 

INTEGER, CONTROL  PARAMETER  THE  NUMBER  OF  INTERVALS  * 
IN  THE  SENSITIVITY  ANALYSIS.  x 


SATELLITE  PHYSICAL  PARAMETERS 


REAL,  MOTOR  SPECIFIC  IMPULSE  IN  SECONDS.  AGRAWAL 
(REF.  2 :  PP .  163-166  )  GIVES  SAMPLE  VALUES. 

.REAL,  INITIAL  MASS  OF  THE  SATELLITE  IN  KILOGRAMS. 
(INCLUDING  THE  FUEL  MASS) 

REAL,  INITIAL  MASS  OF  FUEL  IN  KILOGRAMS. 

REAL,  DRAG  COEFFICIENT  OF  THE  SATELLITE. 
DIMENSIONLESS.  NORMALLY  VARIES  FROM  ONE  TO 
AROUND  THREE,  DEPENDING  ON  VELOCITY,  DENSITY, 
COMPOSITION  OF  ATMOSPHERE  AND  F^YSICAL  DESIGN. 

IS  MOST  OFTEN  SET  AT  2.20  FOR  ALTITUDES  BETWEEN 
200  AND  400  KM.  FOR  ALTITUDES  OF  1000  KM  CD  IS 
SET  AT  JUST  ABOVE  3.0.  JACCHIA  (REF.  3 :P  2). 
.REAL,  SATELLITE  CROSSESTIONAL  AREA  IN  SQUARE 
KILOMETERS. 


ATMOSPHERIC  SPECIFIC  PARAMETERS 


REAL,  RATIO  OF  ATMOSPHERIC  TO  EARTH  ANGULAR  X 
ROTATION  RATES.  A  CONSTANT  IN  EQUATION  11  OF  X 
NSWC  TR  83-243.  UNITLESS.  THE  ACTUAL  VALUE  WILL  * 
DEPEND  ON  LATTITUDE  AND  ALTITUDE.  RANGE  .8  TO  1.0.x 
REAL,  DECIMETRIC  SOLAR  FLUX.  JACCHIA  (REF  4  S  5).  * 
THIS  VALUE  IS  USED  BY  THE  INDUSTRY  TO  REPRESENT  * 
THE  XUV  FLUX  IN  THE  THERMOSPHERE.  THE  XUV  FLUX  * 
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X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 
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IS  THE  FACTOR  CREATING  THE  SUN  FACING  DIURNAL 
DENSITY  AND  TEMPERATURE  BULGE  IN  THE  THERMOSPHERE. 
IT  IS  DIRECTLY  RELATED  TO  SOLAR  ACTIVITY.  IT 
VARIES  WITH  THE  27  DAY  SLOAR  ROTATION  AND  11  YEAR 
SUNSPOT  ACTIVITY  CYCLE.  ORR  (REF.  6:APPENDIX  B) 
INDICATES  RANGES  OF  73.3  TO  242.9.  XUV  FLUX  AND 
THE  GEOMAGNETIC  HEATING  EFFECTS  CAN  CAUSE 
DENSITY  TO  CHANGE  3Y  THREE  ORDERS  OF  MAGNITUDE  AT 


600  KILOMETERS.  THE  10.7  CM  FLUX  IS  DECTABLE  FROM  X 


FBAR 


THE  GROUND  AND  THEORETICALLY  PARELLELS  THE  XUV 
FLUX.  THE  XUV  FLUX  IS  NOT  MEASUREABLE  FROM  THE  x 
GROUND  DUE  TO  NEAR  COMPLETE  UPPER  ATMOSPHERE  x 

ABSORPTION.  x 

REAL,  THE  AVERAGE  10.7  CM  FLUX.  THE  EFFECTS  OF  THEx 
XUV  HEATING  DO  NOT  RETURN  TO  THE  ORIGINAL  UNHEATEDx 


AKP 


TVE 


CONDITION  IN  AN  EARTH  ROTATION.  A  HIGH  HEATING 
CONDITION  WILL  HAVE  SOME  LONGER  TERM  EFFECTS. 
THESE  LONGER  TERM  EFFECTS  ARE  MODELED  HERE  BY 
TAKING  A  MOVING  90  DAY  10.7  CM  AVERAGE. 

FBAR  RANGES  FROM  A  LOW  OF  73  TO  A  HIGH  OF  230. 
REAL,  GEOMAGNETIC  ACTIVITY  INDEX,  JACCHIA  (REF.  4 
AND  5).  GEOMAGNETIC  ACTIVITY  ALSO  CAUSES  HEATING 
IN  THE  THERMOSPHERE.  THE  EXACT  MECHANISM  IS  NOT 
UNDERSTOOD.  IDEALLY  QUEIT  MAGNETIC  CONDITIONS 
CORRESPOND  TO  AN  INDEX  OF  0.  MAXIMUM  ACTIVITY  FOR 
WHICH  SIGNIFICANT  DATA  EXISTS  ARE  RELATED  TO  AN 
INDEX  OF  BETWEEN  6  AND  7 .  NORMAL  ACTIVITY 
INDICATES  AN  INDEX  OF  BETWEEN  1  AND  2. 

REAL,  TIME  OF  VERNAL  EQUINOX  PASSAGE.  USED  TO 
TRANSFORM  TIME  FROM  INERTIAL  TO  EARTH  FIXED 


REFERENCE  FRAMES.  FOUND  IN  THE  ASTRONOMICAL  ALMANAC 


FOR  THE  YEAR  OF  CHOICE.  (MODIFIED  JULIAN  DAYS) 
A  DESCRIPTION  IS  IN  TAFF  (REF.  7:PP.  103-104). 


NAVSPASUR  ORBITAL  ELEMENTS 


THE  NAVSPASUR  ORBITAL  ELEMENT  SET  IS  THE  BROUWER  MEAN 
ELEMENT  SET.  BROUWER  (REF.  1)  IS  THE  BEST  DESCRIPTION  FOR 
THIS  ELEMENT  SET. 


UJD 


ETU 


REAL,  TIME  OF  EPOCH  IN  MODIFIED  JULIAN  DAYS. 

THIS  IS  THE  TO  IN  THE  EQUATION  FOR  THE  MEAN  ANOMALY 
M.  (  M  =  MO  +  N  x  (T  -  TO)  ).  TAFF  (REF.  7).  x 
REAL,  INITIAL  MEAN  ANOMALY  MO.  MEASURED  IN  DEGREES* 


GO. 


HO 


ES. 


HAS  A  RANGE  OF  0  TO  360. 

REAL,  MEAN  ARGUMENT  OF  PERIGEE.  RANGES  FROM  0  TO 
"360  DEGREES.  NOT  DEFINED  FOR  CIRCULAR  ORBITS. 
REAL,  MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE. 
RANGES  FROM  0  TO  360  DEGREES.  NOT  DEFINED  FOR 
EQUITORIAL  ORBITS 

REAL,  MEAN  ECCENTRIY.  FOR  THIS  THEORY  RANGES  FROM 
'0.01  TO  0.1.  THESE  ARE  LOW  ECCENTRICITY  ORBITS. 


XI 


SEAL,  MEAN  INCLINATION.  RANGES  FROM  0  TO  180  x 
DEGREES.  VALUES  OF  0  AND  130  ARE  IN  EQUITORIAL  x 
ORBITS.  INCLINATION  OF  90  IS  A  POLAR  ORBIT.  RANGES* 


RND. 


ESD. 
AO _ 


OF  0  TO  90  ARE  IN  PROGRADE  ORBITS.  INCLINATIONS 
90  TO  180  ARE  IN  RETROGRADE  ORBITS.  FOR  AN 
INCLINATION  OF  0  OR  180  HO  IS  NOT  DEFINED.  THERE 
IS  NO  ASCENDING  NODE  IN  EQUITORIAL  ORBITS. 

REAL,  RATE  OF  CHANGE  OF  MEAN  MOTION.  MEASURED  IN 
'REVOLUTIONS  PER  (DAY  SQUARED) .VALUES  OF  .00001 
AND  ABOVE  ARE  CONSIDERED  SIGNIFICANT  WHEN 
ATMOSPHERE  INDUCED.  REFERENCE  (TAFF)  P137 
REAL,  ECCENTRICITY  DECAY  RATE  MEASURED  IN 
"INVERSE  SECONDS. 

REAL,  KAULA  SEMIMAJOR  AXIS.  MEASURED  IN 
MEAN  EQUITORIAL  EARTH  RADII  (6378.165  KM).  RANGES 
FROM  1.02594  TO  1.16847  TO  MAINTAIN  AN  ALTITUDE  OFx 
100  TO  1000  KM  FOR  RADIUS  AT  PERIGEE.  THIS  MUST  BE* 


CAREFULLY  SCALED  TO  ECCENTRICITY  TO  PREVENT 
CRASHING  THE  SATELLITE.  PERIGEE  ALTITUDES  ABOVE 
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,**  A 


a 


fi 


39 


3 


a 

.v 


-Vfc.  *"• 

■■  wv 


■  -  ^AAi.-Air, 


*  600  KM  ARE  ESCENTIALLY  AVOVE  THE  ATMOSPHERE.  * 

*  ADOT _ REAL,  SEMIMAJOR  AXIS  DECAY  RATE.  * 

X  x 

X  x 

X  * 

X  x 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

DIMENSION  SNSARC  20 ) ,  SNSMIN(20>,  IRVC10),  DAC10),  SNSDEL ( 20 ) 

x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  THE  SET  OADJ  IS  COMMON  TO  PLEP,  K02AI,  THRST  AND  SNALYS.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

COMMON  /  OADJ/  NOA,  IRV,  DA.  SPI,  WT,  PSENS 
X 

.  XXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X  THE  SET  XXXX  IS  COMMON  TO  PLEP  AND  PERIOD.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

COMMON  /  XXXX  /  IET 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  THE  SET  SOLAR  IS  COMMON  TO  PLEP  AND  DRAG.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

COMMON  /  SOLAR  /  F107  ,  FBAR  ,  AKP  ,  TVE 
INTEGER  PSIZE,  PSENS,  PSYNS 
REAL  MINSN 
x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  INITIALIZING  CONSTANTS.  * 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

DATA  R2  /  541 . 15E-06  /  , DEG  /57 . 295779513D0  /  ,  ROE  /6378.165D0/ 
x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  INPUTTING  CONTROL  PARAMETERS.  x 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 

READ  x,  ITER1 ,  ITER2,  NOA,  PSENS,  IET,  PSfZE 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  INPUTTING  SATELITE  SPECIFIC  PARAMETERS.  x 

KXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 

READ  X,  CD,  A,  D,  SPI,  W,  WT 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  INPUTTING  EARTH  AND  ATMOSPHERIC  PARAMETERS.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

READ  x  ,  FI07  ,  FBAR  ,  AKP  ,  TVE 
x 

xxxsrxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  INITIALIZING  BROUWER  MEAN  ORBITAL  ELEMENTS  FROM  THE  NAVSPASUR  x 
X  FIVE  CARD  FORMAT  ORBITAL  DATA  SET.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

READ  5  ,  UJD  ,  ETU  ,  HO  ,  GO  ,  ES  ,  XI 
5  FORMAT  (  /8X,F14.8,5C1X,F8.4)) 

READ  10  ,  RND  ,  ESD 
10  FORMAT  (  20X,  F11.9,  19X,E12.5  ) 

READ  15  ,  AO  ,  ADOT 
15  FORMAT  (  /8X,F11.5,1X,E14.7) 
x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  INPUTTING  THE  REVOLUTION  NUMBER  AND  CORRESPONDING  CHANGE  IN  x 

X  SEMIMAJOR  AXIS  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

IF  (  NOA  .EQ.  0)  GO  TO  30 
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DO  20  K  =  1,  NOA 

READ  *,  IRV(K)  ,  DACIO 
20  CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXXXX 
X  THIS  NEXT  SERIES  OF  CODE  SETS  THE  INITIAL  VALUES  FOR  THE  SCANNING  x 
X  ARRAY  MATRIX  SO  AS  TO  ALLOW  THE  RECOVERY  OF  THE  ORIGINAL  INPUT  x 
x  ELEMENTS  IF  SENSITIVITY  MODE  IS  NOT  SELECTED.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

30  DO  90  ICN  =  1,20 

SNSDEL ( ICN )  =  1.0 
SNSMINC ICN)  =0.0 
SNSARCICN)  =  0.0 
90  CONTINUE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

*  THE  NEXT  SERIES  OF  WRITE  STATEMENTS  OUTPUTS  THE  ORIGINAL  INPUT 

X  PARAMETERS.  UNITS  ARE  SUPPLIED  FOR  CLARIFICATION.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

WRITEC6 , 98 ) 

WRITE(6 , 99 ) 

WRITEC6 ,103) 

WRITE(6,99) 

WRITE! 6 , 1 04 )  D,  F107,  FBAR,  AKP,  TVE 
WRITE(6,I05)  CD,  A,  SPI,  W,  WT 
WRITEC6 , 106 )  UJD,  ETU,  HO,  GO,  XI 
WRITE(6 , 107 )  RND,  ES,  ESD,  AO,  ADOT 
WRITEC6, 103)  ITER1 ,  NOA 
IFCNOA  .GT.  0)  THEN 
WRITEC6, 109) 

DO  27  NN  =  1,  NOA 

WRITEC6, 110)  IRV(NN) ,  DA(NN) 

27  CONTINUE 

END  IF 
WRITE(6 , 99 ) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

*  THIS  NEXT  CODE  SETS  THE  INTIAL  INPUT  PARAMETERS  TO  SOME  RETAINER* 

*  VARIABLES  THAT  ALLOW  RETURNING  TO  THE  ORIGINAL  VALUES  IF  THE  * 

X  SENSITIVITY  MODE  IS  SELECTED.  THIS  PREVENTS  HAVING  TO  REREAD  ALL  * 
x  THE  INPUT  PARAMETERS  BACK  IN  FOR  EACH  RUN  OF  THE  SENSITIVITY  SCAN.* 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXX 

VI  =CD 
V2  =  A 
V3  =  D 
V4=SPI 
V5=W 
V6=WT 
V7  =  F1 07 
V3=FBAR 
V9=AKP 
VI 0= TVE 
V11=UJD 
VI 2= ETU 
V13=H0 
V 1 4  =  G  0 
V15=ES 
V16  =XI 
VI 7  =RND 
V13=ESD 
V 1 9  = AO 
V20  =ADOT 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  THIS  IS  THE  ITERATIVE  SCANNING  ROUTINE.  ARGUMENTS  ARE  INCREMENTED  * 
x  FRACTIONALLY  THROUGH  THEIR  SPECIFIED  RANGE  AS  SPECIFIED  IN  THE  * 
X  INPUT  TO  THE  SUBROUTINE  SNALYS  (PSIZE  ITERATIONS  OF  PLEP).  * 

x  IF  SENSITIVITY  MODE  IS  NOT  SELECTED  PSIZE=1.  * 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


*  * 


'wumtwwiwuiwHininmMMW»M 


IFCPSENS  . EQ .  0)  THEN 
PSIZE=0 
END  IF 

DO  40  PSYNS  3  0,  PSIZE 

x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  THIS  NEXT  SET  OF  CODE  IS  ENVOKED  WHEN  PSENS  IS  EQUAL  TO  ONE.  * 

X  IT  SETS  UP  FOR  AND  CALLS  THE  SENSITIVITY  ARGUMENT  ARRAY  SNALYS  FORX 
X  SUBSEQUENT  SCANNING.  (PSENS  3  1  SELECTS  SENSITIVITY  MODE  * 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

IF  (PSENS  .EQ.  1)  THEN 
IET30 
x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  THE  SUBROUTINE  SNALYS  INTIALIZES  THE  SCANNED  INPUT  PARAMETERS  TO  X 
X  BE  INCLUDED  IN  THE  STUDY.  IT  ALSO  PRINTS  THE  SENSITIVITY  MATRIX.  X 
' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

‘  X 

IF(PSYN5  .EQ.  0)  THEN 

CALL  SNALYS(  PSIZE,  SNSAR,  SNSMIN,  SNSDEL,  SUMCK,  DELSN, 
1  MINSN) 

ELSE 

CD  3  VI 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  RECOVERY  OF  THE  ORIGINAL  INPUT  PARAMETERS  AFTER  THE  FIRST  RUN  OF  x 
X  A  SENSITIVITY  ANALYSIS.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

A  *  V2 
D  3  V3 
SPI3  V4 
W  3  V5 
WT  3  V6 
F107 3  V7 
FBAR=  V8 
AKP  =  V9 
TVE=V10 
UJD=V1I 
ETU=V12 
HO  =V13 
GO  3V14 
ES  =V15 
XI  =V16 
RND=V17 
ESD=V18 
AO  =V19 
ADOT =V20 
END  IF 
END  IF 

ASYNS=REAL (PSYNS) 

ASIZE-REAL( PSIZE) 

IF(PSENS  .EQ.  1)  THEN 

I F( SUMCK  .EQ.  1.0)  THEN 

SCNFRC  =  ASYNSXDELSN  +  MINSN 

ELSE 

SCNFRC  3  ASYNS/ASIZE 
END  IF 
END  IF 

XXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X  THIS  NEXT  CODE  CREATES  THE  CURRENT  VALUE  FOR  THE  SCANNED  INPUTS.  x 
X  IF  SENSITIVITY  IS  SELECTED  THE  PARAMETER  WILL  SCAN  FROM  MIN  TO  MAXX 
X  WITH  A  GRANULARITY  OF  1/PSIZE  PER  SCAN  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

CD  3  SNSAR(1)X(SNSMIN(1)+ASYNSXSNSDEL(1))-CDX(SNSAR(1)-1.0) 

A  3  SNSAR(2)X(SNSMIN(2)+ASYNSXSNSDEL(2) )-A  x( SNSAR( 2) -1 . 0 ) 

D  3  SNSAR(3)X(SNSMIN(3)+ASYNSXSNSDEL(3))-D  x(SNSAR(3)-l .0) 

SPI3  SNSAR(4)X(SNSMIN(4)+ASYNSXSNSDEL(4))-SPIX(SNSAR(4)-1 .0) 
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W  *  SNSAR(5)»(SNSMIN(5>+ASYNSXSNSDEL(5))-W  X(SNSAR(5)-1 . 0) 

WT  *  SNSAR(6)X(SNSMIN(6)+ASYNSXSNSDEL(6) )-WTx( SNSARC  6  )-l . 0  ) 

F107  =  SNSARC7)XCSNSMIN(7)+ASYNSXSNSDELC7))-F107x( SNSARC 7  )-l .0) 
FBAR*  SNSAR(S)x(SNSMIN(8 )+ASYNS*SNSDEL (3))-FBARX(SNSAR(8)-l  .0) 

AKP  =  SNSARC  9 )X( SNSMINC 9 ) +ASYNSXSNSDEL (9))-AKPx(SNSAR(9)-l.Q) 
TVEsSNSAR(10)x(SNSMIN(10)+ASYNS*SNSDELC10))-TVEx(SNSAR(10)-l  .0) 
UJD3 SNSARC ll)x(SNSMIN(ll)+ASYNS*SNSDEL(ll))-UJD*(SNSAR(ll)-l  .0) 
ETUsSNSAR(12)*(SNSMIN( 12 )+ASYHSXSNSDEL (12))-ETUX($NSAR(12)-1 .0) 

HO  =SNSAR( 13 )x( SNSMINC 13)+ASYNSXSNSDELC 13) )-HQX( SNSARC 1 3 ) - 1 .0) 

GO  =SNSAR( 14 )x( SNSMINC 14 )+ASYNSXSNSDELC 14) )-GQx( SNSARC 14) -1 .0) 

ES  =SNSAR(15)X( SHSMINC 15)+ASYNS*SNSDEL  (15))-ESx( SNSARC 15  )-l  .  0 ) 

XI  *SNSARC 16 )X( SNSMINC 16 )+ASYNSXSNSDEL C 16 ))-Xlx( SNSARC 165-1.0) 
RND=SNSAR(17)X(SNSMIN(17)+ASYNS*SN5DEL(17))-RNDX( SNSARC 1 7 ) — 1 .0) 
ESD=SNSARC 18 )x( SNSMINC 18 )+A5YNSXSNSDEL (18))-ESD*(SNSAR(18)-1 .0) 

AO  =SNSARC 19)x( SNSMINC 19)+ASYNSXSNSDELC 19) )-AOx( SNSARC 19 )-l .0) 
ADOT*SNSARC 20) XC SNSMINC 20 )+ASYNSXSNSDELC 20 ))-ADOTXC SNSARC 20 )-l .0) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  DEL  IS  THE  CONSTANT  MULTIPLIER  FOR  THE  INTEGRAL  IN  EQUATION  33  OF  x 
X  N5WC  TR  83-243  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

DEL  =  CDXAX3. 14159D0/SPI  . 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  IF  THE  SENSITIVITY  ANALYSIS  MODE  IS  NOT  SELECTED,  THEN  THE  SYSTEM  x 
X  PRINTS  STANDARD  MISSION  LIFE  TIME  TABULAR  DATA.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

IF  CPSENS  .EQ.  0)  THEN 
BCOF  *  CD  x  A 
WRITE  (6,98) 

WRITE  C 6 , 99) 

WRITE  (6,100) 

WRITE  (6,99) 

WRITE  C 6 , 1 0 1 )  BCOF,  SPI,  W 
WRITE  (6,102)  F107 ,  FBAR,  TVE 

END  IF 

x 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXX 

X  CONVERTING  DEGREES  TO  RADIANS  » 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXX 

X 

XI  a  XI  /  DEG 
GO  =  GO  /  DEG 
HO  a  HO  /  DEG 
FLO  *  ETU  /  DEG 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  THIS  NEXT  CODE  SEQUENCE  CONDITIONS  THE  ELEMENTS  AS  PER  NSWC  TR  82-387 
X  PAGE  7.  IT  CONVERTS  THE  KAULA  SEMIMAJOR  AXIS  AO  IN  EARTH  RADII  TO  x 
X  BROUWER  MEAN  SEMIMAJOR  AXIS  VIA  EQNS  2.1  AND  2.2.  THE  SEMIMAJOR  AXISX 
X  DECAY  RATE  ADOT  IS  CONVERTED  IN  A  LIKE  MANNER  VIA  2.3,  2.4,  AND  2.5  x 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

FN  =  1 .  /  (  1 .  -  ES  x  ES  )XX1 .5D0 
TA  =  SIN  (  XI  ) 

TB  =  TA  x  TA 

TE  a  1.  -  (3.XTB)  /  2. 

TAD  «  ((  3.XR2)  /  C2.XA0XA0))XTEXFN 
AK  a  AO 

AO  =  CA0XCC1.  +  2.XTAD)  /  (1.  -  TAD) )xxo . 66666667D0 ) xROE 
XDOT  =  TADXC3.XESDXCES/C1.-ESXES))  -  2 . XCADOT/AK) ) 

ADT  »  (ADOTZAK)XAO  +  2 . XAKXROEXXDOTXC ( 1 .  ♦  2.xTAD)x 
1  C(l.  -  TAD) xx  5))xx(-0. 3333333333 DO) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  TIME  IS  NOW  CONVERTED  FROM  SECONDS  TO  DAYS  * 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


ADT  «  ADT  /  86400. DO 
ESD  a  ESD  /  86400. DO 


RND  *  RNDX2.D0X3.14159D0X(1./(86400.D0)XX2) 
x 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X  THE  SUBROUTINE  KOZAI  IS  THE  OPERATING  HORKHORSE  FOR  THE  PROPELLENT  x 
X  LONGEVITY  COMPUTATIONS.  IT  ALSO  DECREMENTS  PROPELLENT  AS  IT  IS  USEDX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

CALL  KOZAK A DT ,AO,ES,XI,FL0,GO,HO,UJD,W,DEL, ITER1 , ITER2, 

1  SCNFRC,  D) 

x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  FORMAT  STATEMENTS  USED  IN  THE  PLEP  INITIALIZATION  PROGRAM.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

98  FORMAT  (  1H1  ) 

99  F0RMAK//5X,  ' - ',//) 

100  F0RMAT(//15X, 'PROPELLANT  LIFE  ESTIMATOR*,//) 

101  F0RMATC//5X, ' CDA2 • , E12 . 5,  *  SQ  KM*/5X, 'SPECIFIC  IMPULSE2 ' , 

XF7 . 2, *  SEC'/SX, 'INITIAL  WEIGHT  OF  PROPELLANT2', 

XF7 . 2, '  KG') 

102  FORMAT (//5X, 'F1072',F8 .2/5X, * FBAR= ' , F8 . 2/5X, 

X  '  TIME  OF  VERNAL  EQUINOX2 ', E22 . 14) 

103  FORMAT ( 1SX,  'INITIAL  INPUT  CONDITIONS') 

104  F0RMAT(///1X, 'ATMOSPHERIC  INPUT  PARAMETERS »//5X, ' D  (RELATIVE  ATMOS 
XPHERIC  ROTATION)  =  ' , F10 . 4/5X, ' F107  (DECIMETRIC  SOLAR  FLUX)  = 

XF7 . 2/5X, ' FBAR  (90  DAY  AVG.  F107)  =  ' , F7 . 2/5X, ' AKP  (GEOMAGNETIC  IND 
XEX)  2  ' ,  F6 .3/5X, 'TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  2  ',F16.9, 

X*  MODIFIED  JULIAN  DAYS'///) 

105  FORMATdX, 'SATELLITE  PHYSICAL  PARAMETERS '//5X, 'CD  (DRAG  COEFFICIEN 
XT)  2  ' , F5 . 2 , '  UNITLESS '/5X, 'A  (CROSSECTIONAL  AREA)  2  *,F11.9,*  SQU 
XARE  KIL0METERS'/5X, 'SPI  (MOTOR  SPECIFIC  IMPULSE)  2  *,F8.2,*  SECOND 
XS’/5X,'W  (INITIAL  FUEL  MASS)  2  *,F9.2,*  KILOGRAMS'/5X, *WT  (INITIAL 
X  SATELLITE  MASS)  2  ',F11.2,'  KILOGRAMS*///) 

106  FORMAT (  IX, ' NAVSPASUR  ORBITAL  ELEMENTS '//5X, 'UJD  (EPOCH)  2  ',F14.8 
X,’  MEAN  JULIAN  DAYS'/5X, 'ETU  (MEAN  ANOMALY J  2  ',F8.4,*  DEGREES V5X 
X, ' HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  2  ',F8.4,*  DEGRE 
XES'/5X, 'GO  (MEAN  ARGUMENT  OF  PERIGEE)  2  ',F8.4,'  DEGREESV5X, 

X'XI  (MEAN  INCLINATION)  =  ',F8.4,'  DEGREES') 

107  F0RMAT(5X, 'RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  2  *,F11.9,'  REVOLUT 
XTIONS  PER  DAY  SQUARED'/ 

X5X, ' ES  (MEAN  ECCENTRICITY)  2  ',F8.4,'  UNITLESS'/5X, ' ESD  (ECCENTRIC 
XITY  DECAY  RATE)  2  ',FI6.9,'  1/SECONDS '/5X, 'AO  (KAULA  SEMI-MAJOR  AX 
XIS)  2  ' , Fll . 5, '  MEAN  EARTH  RADII ’/5X, *ADOT  (SEMI-MAJOR  AXIS  DECAY 
XRATE)  2  ' , F16 . 9 , '  •///) 

108  FORMATdX, 'CONTROL  INPUT  PARAMETERS '//5X, 'MAXIMUM  REVOLUTIONS  PER 
XITERATION  17 ,/5X, 'NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJO 
XR  AXIS  ' , 12/  ) 

109  F0RMAT(//5X, 'SEMI -MAJOR  AXIS  INCREASE  SCHEDULE'//5X, 'REV  ADJUST', 
X5X, 'CHANGE  IN  KILOMETERS'/) 

110  FORMAT (5X,I7,7X,F6.2) 


x 

40  CONTINUE 
:  OP 
END 
x 

X 

X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  xx  SUBROUTINE  SNALYS  xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

X 

SUBROUTINE  SNALYS  (PSIZE, SNSAR,  SNSMIN,  SNSDEL,  SUMCK,  DELSN, 
1MINSN ) 


DIMENSION  SNSAR(20),  $NSMIN(20) ,  SNSMAX(20),  SNSDELC20) 

CHARACTERX4  SNSCAR 
INTEGER  N,  PSI2E 

REAL  SNSAR,  SNSMIM,  SNSMAX,  SNSDEL,  SUMCK,  DELSN,  MINSN 
WRITE  (6,98) 

WRITE(6 , 120 ) 

WRITE(6, 100) 

WRITE(6 ,120) 

WRITEC6, 200 )  PSIZE 
ASIZE=REAL(PSIZE) 

READC 5,102) 

SUMCK  =0.0 
DO  20  N  =  1,  20 

READ  101,  SNSCAR,  SNSAR(N) ,  SNSMIN(N),  SNSMAX(N) 

SNSDEL(N)  =  (SNSMAX(N)-SNSMIN(N))/ASIZE 
IF  (SNSAR(N)  .EQ.  1.0)  THEN 
SUMCK  =  SUMCK  +1.0 
DELSN  =  SNSDEL(N) 

MINSN  =  SNSMIN(N) 

WRITE(6 , 300 )  SNSCAR,  SNSMIN(N),  SNSMAX(N), 

1  SNSDEL(N) 

END  IF 
20  CONTINUE 

IF  (SUMCK  .GT.  1.0)  THEN 
DELSN  =0.0 
MINSN  =0.0 
END  IF 

HRITE(6 , 120 ) 

WRITE(6 , 400 ) 

98  FORMAT  (  1H1  ) 

120  F0RMAT(//5X,  * - ',//) 

100  FORMAT(//17X, 'SENSITIVITY  ANALYSIS  FOR' ,/16X, 'LOW  EARTH  ORBIT  ', 

X  'SATELLITES', /19X, 'PROPELLENT  LONGEVITY',//) 

200  FORMAT  (5X,  •  TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO', 

X  /5X, '  MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN.  ',/5X, 

X  •  THE  GRANULARITY  IS  1/ • , 14, ' . */5X,  *  ALL  OTHER  INPUT  PARAMETERS 
X  ARE  RETURNED  TO  THEIR', /5X, '  INITIAL  VALUES  AT  EACH  INCREMENT.' 

X  ///'INPUT  PARAMETER' ,8X, 'MINIMUM* , 13X, 'MAXIMUM* , 13X, • INCREMENT'/) 

102  FORMAT  (IX) 

101  FORMAT  (  A4 ,  IX,  F3.1,  F16.9,  4X,  F16.9) 

300  FORMAT  (5X,  A4,  6X,  3(3X,  F16.9)) 

400  FORMAT(// ,  10X,  '  SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS'// 
15X, 'REV  NUMBER', 8X, 'TIME( DAYS) • , 15X, 'FRACTION  OF  SCAN') 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

xx  SUBROUTINE  KOZAI  xx 

XX 
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X 

SUBROUTINE  KOZAI  IS  THE  SYSTEM  WORKHOUSE.  IT  CALLS  BRAUER  TO  CON-x 
VERT  THE  BROUWER  MEAN  ORBITAL  ELEMENTS  OF  THE  NAVSPASUR  INPUT  TO  x 
BROUWER  OSCULATING  ELEMENTS.  THESE  OSCULATING  ELEMENTS  ARE  CON-  X 
VERTED  TO  KOZAI-lIKE  MEAN  ELEMENTS  IN  THE  SUBROUTINE  MEAN.  THE  x 
SUBROUTINE  PERIOD  COMPUTES  THE  INITIAL  ANOMALISTIC  PERIOD.  THE  * 
SUBROUTINE  GEOP  COMPUTES  THE  EFFECTS  OF  THE  OBLATE  GEOID  ON  THE  x 
ORBITAL  ELEMENTS  (THROUGH  THE  J4  TERM).  SUBROUTINE  THRST  COMPUTES  x 
THE  EFFECTS  ON  THE  ELEMENTS  DUE  TO  PERIGEE  BURNS  TO  CHANGE  x 
THE  SEMIMAJOR  AXIS  DURING  IRV(K)  BY  THE  AMOUNT  DACK).  FINALLY,  X 
SUBROUTINE  DRAG  COMPUTES  THE  THE  AMOUNT  OF  FUEL  NEEDED  TO  OVERCOME  x 
THE  DRAG  EFFECTS  3Y  CONTINUOUS  INTRACK  MICROTHRUSTI NG .  X 
KOZAI  KEEPS  TRACK  OF  THE  FUEL  USED  AND  ORBITAL  ELEMENT  CHANGES.  x 
WHEN  THE  FUEL  IS  EXHAUSTED  KOZAI  RETURNS  SYSTEM  CONTROL  TO  PLEP.  x 
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SUBROUTINE  KOZAI ( ADT, AO , ES, XI , FLO, GO , HO, UJD, W, DEL , ITER1 , ITER2, 

1  SCNFRC,  D) 

DIMENSION  XOSCC6)  ,  XM(6),  DXDTC6)  ,  XDDC6),  IRV(IO),  DA(IO) 
DIMENSION  DXMC6) 

COMMON  /  OADJ/  NOA,  IRV,  DA,  SPI,  WT,  PSENS 

DATA  BO, B2, B3, B4, B5, B6  /  . 3986 03254E+06 , - . 1755528999E+11 , 

1  . 26 386 6 47 738 E+l 2, . 1063073996E+16, .80560S022E+18, 

1  . 7292115856E-04/ 

DATA  DEG/57 . 295779513DQ/,CN2/0 . 0/, A1 , El , RN1/Q . 0, 0 . 0, 0 . 0/, DT/O . 0/ 
DATA  PI  /  3 . 14159D0  / 

DATA  NPLOT  /  0  / 

IRC  =  0 
PC0UNT=0 
KK  =  0 

PI2  =  2.D0X  PI 

CALL  BRAUER  ( BO , B2. B3 , B4 , 85, DT , AO , ES,XI , FLO, GO ,H0 ,CN2, A,CE, Cl , CL , 
1  G,H, A1 , El , RN1  ) 


IF  (PSENS. EQ.  0)  THEN 
WRITEC6 ,100)  ITER1 

100  F0RMAT(//5X, 'MAXIMUM  NUMBER  OF  ORBITAL  REVOLUTIONS* ', II 0 ) 
WRITEC6 , 101 )  ITER2 

101  F0RMAT(//5X,  'DATA  PRINTED  EVERY  M3,'  ORBITAL  REVOLUTIONS') 
SCNFRC  *  0.0 

END  IF 

CALL  PERIOD  (  AO,  ES,  XI,  PA  ) 

XOSC(l)  =  A 
X0SCC2)  =  CE 
XOSC( 3 )  =  Cl 
X0SCC4)  *  G 
XOSC( 5  )  =  H 
XQSCC6)  =  CL 
CALL  MEAN  (XOSC.XM) 

ILINE  =  0 
CIW  =  XM(3)XDEG 
GW  =  XM(4)XDEG 
HW  s  XM( 5 ) XDEG 
CLW  =  XM( 6 ) XDEG 
IF  (PSENS  .EQ.  0)  THEN 
WRITEC6 ,102) 

WRITE(6 ,103)  (XM(K),K*1,2),CIW,GW,HW,CLW 
END  IF 

103  F0RMATC/5X, ’SEMIMAJOR  AXIS** , E22 . 14, •  KM’/ 

X  5X, 1  ECCENTRICITY2 ' , E22 . 14/ 

X  5X, ’INCLINATION*’ ,E22. 14,  '  DEG’/ 

X  5X, ’ARGUMENT  OF  PERIGEE* ' , E22 . 14,  '  DEG’/ 

X  5X, 'ASCENDING  NODE*' , E22 . 14,  '  DEG’/ 

X  SX, 'MEAN  ANOMALY*', E22. 14, '  DEG') 
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102  F0RMAT(//5X, 'KOZAI  MEAN  ELEMENT  SET*) 

DELT  =  PA 

CALL  THRST  (  XM,  DXM,  IRC,  KK,  PA,  DWOA) 

DO  10  1=1,  ITER1 

T  =  (I-1)X  DELT  /  86400. 

TT  =  UJD  +  T 
CALL  GEOP  (XM.DXDT) 

DO  20  JJ  =  1,6 

XM(JJ)  =  XMCJJ)  +  DXDT ( J J )*DELT  +  DXM(JJ) 
IF  (  JJ  .  LT  .  4  )  GO  TO  20 

XM(  JJ  )  =  AMOD  (  XMCJJ),  PI2  ) 

CONTINUE 
IRC  =  0 

DO  30  II  =  1,  NOA 

IF  (  I.EQ.  IRV(II))  IRC  =  1 
IF  (  IRC  .EQ.  1)  KK  =  II 
CONTINUE 

CALL  THRSTCXM,  DXM,  IRC,  KK,  PA,  DWOA) 


20 


30 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X  CALCULATING  THE  DELI  IN  NSWC  TR  83-243  EQN  11.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

RP=XM(l)x(l. O-XM! 2) ) 

VP=SQRT(CB0/XMC1))X((1 . 0+XMC2) )/( 1 . 0-XM(2) ) ) ) 

DEL1=(1 . 0-(RP/VP)XDXB6XCOS(XM(3)))XX2 
DDELT=DEL*DEL1 


CALL  DRAG  (  XM,  DDELT,  TT,  DELMG  ) 

W  =  W  +  DELMG  -  DMO A 
WT  *  WT  +  DELMG  -  DMOA 
IF  (PSENS  .EQ.  0)  THEN 

IF  (PCOUNT.NE. CITER2-1))  THEN 
PCOUNT =PC0UNT+1 
GO  TO  10 

ELSE 

PCQUNT=0 
END  IF 

IF  <  ILINE  .GE.  50  )  ILINE  =  0 
ILINE  =  ILINE  +  1 

IF  (  ILINE  .EQ.  1  )  WRITE  (  6,104) 

IF  (  ILINE  .EQ.  1)  WRITE  (  6,105) 

WRITE  (  6,  106)  I,  TT,  T,  W 
END  IF 

IF  (  W.LE.  0.00)  THEN 

IF  CNPLQT  .EQ.  1)  THEN 
PRINT  X,  T,  SCNFRC 

ELSE 

WRITEC6 ,107)  I,  T,  SCNFRC 
END  IF 
RETURN 
END  IF 
10  CONTINUE 

IF  (NPLOT  .EQ.  1)  THEN 
PRINT  X,  T,  SCNFRC 

ELSE 

WRITEC6 , 107 )  I, T, SCNFRC 
END  IF 

104  FORMAT ( 1H1 ) 

105  F0RMATC/5X, 'REV  NUMBER ', 3X, 'TIME(MJD) • ,8X, 'TIME(DAYS) • ,8X, 
X  'PROPELLANT  REMAINING! KG) ' , // ) 

106  FORMAT (5X,I10,3X,F9.2,8X,F10.2,15X,F8.2) 

107  FORMAT ( 5X, I 1 0 ,  8X, F10 . 2, 15X, F16 . 9 ) 

RETURN 

END 
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SUBROUTINE  BRAUERC BO , B2, B3 , 84 . B5, DT, A2P, E2P, CI2P, CL 02P, G02P , H02P, 
1CN2,A,CE,CI,CL,G,H,ADT,RND,ESD) 

1004  A2P2=A2P**2 

A2P6sA2P2XX2 

CN0*SQRT ( BO/C A2P2XA2P ) ) 

E2P2*E2P**2 
ETA*SQRT ( 1 . -E2P2) 

SINEI=SIN(CI2P) 

THETA*C0S(CI2P) 

TH£TA2aTHETA**2 
THETA4=TH£TA2**2 
THETA6* THETA 4XTHETA2 
CJ2=-B2/( 2 . XB0XA2P2) 

ETA2=ETA**2 

ETA3*ETA2*ETA 

ETA4*ETA2**2 

CJ21P*Cvl2/ETA4 

CJ31P=S3/<B0XA2P2*A2P*E7A4*ETA2) 

CJ41PS(3.*B4)/(8.*B0*A2P4*ETA4*ETA4) 

CJ51P=B5/(B0*A2P4XA2P*ETA4X*2*ETA2) 

FUN 1*3 . XTHETA2-1 . 

FUN2=1 . -5 . XTHETA2 

SINEI2*SINEI**2 

A1=A2P*CJ2*FUN1 

A0*-A1/ETA3 

A2a3 . *A2P*CJ2*SINEI2 

FUN5*1.-11.*THETA2-(40.*THETA4)/FUN2 

FUN6  =  -FUN1-(S.*THETA4)/FUN2 

FUN4=THETA2/SINEI2 

FUN22=FUN2**2 

CJ21P2=CJ21PXX2 

E01P*< (E2P*ETA21*<3 . *CJ21P2*FUN5-10 . * 

1CJ41P*FUN6))/(24.*CJ21P) 

E21P=-2 . xEOlP 

E31P=< (35.*CJ51P*E2P2*6TA2*SINEI)*(FUN2-(16 . *THETA4)/FUN2) )/C96 . * 
1CJ21P) 

E11P*-.75D0xE31P+(( . 25D0*ETA2*$INEI JXCCJ31P+ . 3125D0*CJ51P 
1*(4.+3.*E2P2)» 

2( 1 . -9 . *THETA2-( 24 . *THETA4 )/FUN2> ) )/CJ21P 
CI0=>-(E2P*THETA)/(ETA2*SINEI) 

CI2»CJ21P*THETA*SINEIX1.5D0 
CIl*E2P*CI2x. 66666667  DO 
FUM7=(-.5D0*ETA3*CJ21P)/E2P 

C121P=(ETA3/CJ21P)*( .25D0*CJ21P2*FUN5-.83333333D0*CJ41P*FUN6) 
CL12P*CN0x(l .  +  1 .5D0*CJ21PXETA*FUNl+.09375D0*CJ21P2xETA*(-15.+16  . 
l*ETA+25.*ETA2+(30 .-96 .XETA-90 . *ETA2 ) XTHETA2+( 1 05 . *144 . XETA+25 . 
2*ETA2)XTHETA4)+ . 9375D0*CJ41P*ETA*E2P2*( 3 . -30 . XTHETA2+35 . *THETA4) ) 
CL22P* .5D0XCN0XCN2 

G21P*(1./(24.XCJ21P);*(-3.*CJ21P2*(2.+E2P2-11 . *( 2 . +3 . *E2P2) XTHETA2 
l-40.*(2.+5.*E2P2)*THETA4/FUN2-400 .*E2P2*THETA6/FUN22)+10.xCJ41?* 
2(2. 

3+E2P2-3 .  *( 2 . +3 . *E2P2 )*THETA2-8 . *<  2 . +5 . *E2P2) *THETA4/FUN2-80 . *E2P2* 
4THETA6/FUN22) ) 

G12P=CNQX(-1 .5D0*CJ21P*FUN2+. 09375D0*CJ21P2*(-35.+24.*ETA+25. 
1XETA2+ 

2(90.-192. X ETA-126 . *ETA2 ) *THETA2+( 385 . +360 . XETA+45 . *ETA2)xTHETA4)+ 
3.3125D0XCJ41P»(21 . -9 . XETA2  +  ( -27 0 , +126 . *ETA2 ) *THETA2+( 385 . -189  . 
4*ETA2USTHETA4)) 

H2*l . 5D0XCJ21PXTHETA 
H3=-2.*H2 

HI = . 66666667D0XE2PXH2 

H31P=((35. *CJ51P*E2P2*E2P*THETA)/C144.*CJ21P) )*( . 5D0/SINElx( FUN2-( 
116 . XTHETA4 )/FUN2)+SINEI*( 5 .+C32.XTHETA2)/ FUN2+80 . XTHETA4/FUN22 ) ) 
H11P=- . 25* 

1  H31P+(( .25D0*E2P*THETAJ/(CJ21P*SINEI) )*(CJ31P+. 3125D0XCJ51P* 

2(4.+3.*E2P2)*(l.-9.*THETA2-(24.*THETA4)/FUN2)+l .875D0*CJ51P 
3XSINEI2X 

4( 4 .+3 . *E2P2)*( 3 . +( 16 . XTHETA2 )/FUN2+( 40 . *THETA4)/FUN22)) 


142 


H21P=(E2P2*THETA)/(12.*CJ21P)K(-3.XCJ21P2*U1 .+(30 . *THETA2)/FUN2+ 
1(200 .xTHETA4)/FUN22)+10.xCJ41Px(3 .+(16 . *THETA2)/FUN2+( 40 . XTHETA4 )/ 
2FUN22) ) 

H12P=CN0*THETAX(-3.*CJ21P+.375D0*CJ21P2*(-5.+12.xETA+9.XETA2+ 

1  ( -35 .  - 

236 . XETA-5 . *ETA2)*THETA2)+1 . 25D0xCJ41Px( 5 . -3 . *ETA2)X(3 . -7 . XTHETA2) ) 
AID=CJ51P/CJ21P 
AID2=FUN2-( 16 . XTHETA4)/FUN2 
C1=35./384.XAID*ETA3XE2PXSINEIXAID2 
AID3=THETA2/SINEI 
AID4=TH£TA2*SINEI 
E2P3=E2P2x£2P 

C2=35./1152.XAIDX((-E2PXSINEIX(3 .+2. *E2P2)+E2P3*AID3 >*AID2+ 

12. »E2P3*AID4X(5 . +( 32 . XTHETA2 )/FUN2+(80 . XTHETA4)/FUN22 ) ) 

C3S1 ,-9.xTHETA2-(24. XTHETA4)/ FUM2 
AID5=CJ31P/CJ21P 

C4».25D0XAID5x(-E2PXAID3)+5./64.*AID*(-E2P*AID3*(4.+3.xE2P2)+ 
1E2Px3INEIx(26 .+9.XE2P2) 5XC3-15 . /32 . *AI D*E2P*AI D4X( 4 . +3 . *E2P2) x 
2(3.+(16.XTHETA2)/FUH2+(40.XTHETA4)/FUN22) 
C5=E2P/(1.+£TA3)X(3.-E2P2X(3.-E2P2)) 

C6  s( E2Px( -32 . +3 1 .x(E2°2xE2P2)))/((4.+3.xE2P2)+ETAX(4.+9.xE2P2)) 
C7=.25D0XAID5x3INE:xC5+5./64.XC3XAIDxETA2XSIHEIxC6 
C8=-.25D0XAID5xETA3xSINEI-S./64.xAIDxETA3xSINEIx(4.+9 .XE2P2)XC3 
0510  T  =  DT 

CL2P=CL 12PxDT+CL22PxDTxx2+CL  02 P  +  RHDxDT*DT 
Cl2P=AM0D(Ct2P,6 .283 185307 17 96D0) 

IF(Ct2P)520,530,530 
520  CL2P=CL2P+6 .28 318530717 96 DO 
530  G2P=G12PxDT+G02P 
H2P*H12PXDT+H02P 
SINEG=SIN( G2P ) 

COSING=COS(G2P) 

D1E=SIHEGX(SINEGX(E31PXSINEG+E21P)+E11P)+E01P 
H1P=( ( H31P*SINEG+H21P)»SINEG+H1 IP ) XC0SIHG+H2P 

GPIP=G2P+C12P+ . 5D0X(CL21P+G21P )XSIN(2 . XG2P)+(C1+C2)XC0S( 3 . x  G2P ) 
1+(C4+C7 1 XCOSING 
CL1P=CL2P 
U=CL2P 

100  DELTAU=(U-E2PXSIN(U)-CL2P)/(1.-E2PXC0S(U)) 

U=U-DELTAU 

IF(ABS(DELTAU)-1. £-10)200. 100, 100 
200  U=U-(U-E2PXSIN(U)-C12P)/(1.-E2PXC0$(U)) 

E  =  U 

SINE1P=SIN(E) 

C0SE1P=C0S( E) 

G1P=G2P 

ADIVR=1 . /( 1 . -E2PXC0SE1P) 

SIHF1P=ADIVRXETAXSINE1P 

C0SF1P=ADIVRX(C0SE1P-E2P) 

F1P*ATAN2(SINF1P,C0SF1P) 

IF(ABS (FlP-CL2P)-3. 1415926535898 DO )220, 210, 210 
210  STOP 

220  FUM3  =  ( 1 .+CN2XT)XX. 6666666 7 DO 
C0SFG=C0S(2.X(G1P+F1P)) 

SINFG=S1N(2.*(G1P+F1P)) 

ADIVR2=ADIVR*X2 

ADIVR3=ADIVR**3 

CI-CI2P+CIQ*D1E+CI1*SINF1PxSINFG+(2.*CI1xCQSF1P+CI2)*C0SFG 
FUf.’8  =  FlP-CLlP  +  E2PXSIHFlP 

8018  H=HlP+(2 .XH1HC0SF1P+H2)XSIHFG-H1XSINF1PXC0SFG+H3XFUN8 
KFUN  =  H/6 .28318 530  7 17 96  DO 
FUH9=KFUN 

H=H-FUN9X6 .28 3185307 1796D0 
IF(H)8022, 8023, 3023 

8022  H=H+6 . 2831S5307 1796D0 

8023  A=A2P/FUN3+AO+(A1+A2*COSFG)XADIVR  +  ADTXDT 
AID6=ADIVR2XETA2+ADIVR 

AID7*SIN(2. XG2P+F1P ) 

AID8=SIN(2.XG2P+3.XF1P) 

D1*.25D0xCJ21Px(6.*(5.xTHETA2-1. )xFUN8+(3.-5.xTHETA2)x(3.xsINFG+ 

13. XE2PXAID7  +E2PXAID8  )) 

D2'.25D0XCJ21PX(2.*(3.XTHETA2-1. )x(AID6+l . )xsiNFlP+3 . x( 1 . -THETA2)x 


B  *  K«  s*  M*TDrKJruTyirw*ru>ruxM  iru  *fv  ■ru 


1  ((-AID6+1 .  )xAID7+(AID6+.33333333D0)xAIDS)) 

AID9*C0S(2 . »G2P+F1P) 

AID10*C0SC2.xG2P+3. XF1P) 

03*-ETA2*.5DC*C.I21Px(l.“THETA2)x(3.xAID9+AID10) 

ETA6I*1 ./(ETA3XETA3) 

D4*ETA6Ix(C5  +C3SF1PX(3.+E2P»C0SF1PX(3. +E2P*CQSF1P ) ) ) 
D5*ETA6IX<E2P+COSF1PX(3.+E2PXCOSF1PX(3.+E2PXCOSF1P))) 

06 *  ETA2XCJ2X. 5D0X((3.xTHETA2-l . ) XD4+3 .*(1 .-THETA2) XD5*CQSFG)+D3 
GAL*GPLP+DH-C£2PXETA2)/<1 .+ETA)*D2 

CE*(E2P-1 . )*(  1 . +CN2»DT )»X.66666666666667D0+1.+D1E+D6  +  ESDxDT 
EDL  * . 5D0»£2P»CL21PXSIN  ( 2 . *G2P ) +C3XC0SING+E2PXC1XC0S  ( 3 . XG2P ) - 
1  ETA3*D2 
AID14*SIN(CL2P) 

AID15=C0S(CL2P) 

ESL  *CE  XAID14  +  EDLXAID15 
ECL *CE  XAID15-EDLXAID14 
CE=$?aT(ECL*ECl  +  ESUESL) 

CL*ATAN2(ESL,£CL) 

G=GAL -CL 

G=AM0D(G,6 .2831853071796  DO) 

IF(G)8024, 3025, 8025 
*  G  *  G  ♦  6 .283135307179600. 

5  RETURN 
END 
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SUBROUTINE  POSVEL  ( A, CE, Cl , CL , G, H, BO , R1 , R2, R3, VI , V2, V3, R, VEL ) 
CALL  NWTRPH  (CL , CE, E, SINEP, CQSEP) 

HSIN=SIN(H) 

HCOS*COS(H) 

GSIN-SIN(G) 

GCOS=COS(G) 

CISIN*SIN(CI ) 

cicos*cosccn 

All *HCQS*GCOS-HSINxCICOS*GSIN 

A12*-HCOSXGSIN-HSINXCICOSXGCOS 

A21*HSINXGC0S+HC0SXCIC0SXGSIN 

A22=HC0SXCIC0SXGC0S-HSINXGSIN 

A31=CISINXGSIN 

A32=CISINXGC0S 

FUN* SORT (1 . -CE*CE) 

R1=AX(A11X(C0SEP-CE)+A12X(FUNXSINEP) ) 
R2=AX(A21X(C0SEP-CE)+A22X(FUNXSINEP>) 
R3=AX(A31X(CQSEP-CE)+A32X(FUNXSINEP) ) 

R=AX( 1 . -CEXC3SEP) 

F'JMl*(S3nT(S0»A))/R 

V1=FUN1X( A1 1X( -SINEP )+A12*FUH*COSEP ) 

V2*FUH1 x( A21 *( -SINEP)+A22XFUHXC0SEP) 

V3  =  Fl,tll*(A3I*<-SIME?)+A32*FUNxC0SEP) 

VEL=FUN1XS0RT<1 .-(CEXCOSEP)XXZ; 

RETURN 

END 
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xx  SUBROUTINE  NWTRPH  xx 
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SUBROUTINE  NWTRPH  (QL , B, E2, SE, CE) 

NEWTON  RAPHSON  SOLUTION  TO  KEPLERS  EQUATION  SUBROUTINE 
E1=QL+(BXSIN(QL) V(1 .-BXCOSCQL)) 

2  SE=SIN( El ) 

C£sC0S(  El ) 

E2=E1+(QL+B*SE-E1)/(1.-BXCE) 

IF( ABS( E2-E1 )-  l.E-8  )4,4,3 

3  E1  =  E2 
GO  TO  2 

4  SE=SINCE2) 

CE=C0S(E2) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

XX  FUNCTION  XSPA  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


THE  FOLLOWING  FUNCTIONS  DETERMINE  THE  FIRST  ORDER  SHORT 
PERIODIC  VARIATIONS  AS  A  FUNCTION  OF  THE  MEAN  ELEMENTS. 

FUNCTION  XSPA(XM) 

DIMENSION  XM( 6 ) 

REAL  J(4) 

DATA  J.R/l. ,  0.0010827600,-0. OQQO0255DO,-O. 0Q000156D0, 6  378 .165D0/ 

F  =  XM(6)+(2.xXM(2)-XM(2)XX3./4+5./96 .XXM(2)XX5. )XSIN(XM(6)) 
1+<1.25D0XXM(2)XX2.-11./24.XXM(2)XX4.+17 . / 192 . XXMC2) XX6 . IxSINC  2 . x 
2XM( 

26))+(13./12.xXM(2)xx2.-43./192.xXM(2)xx5.)xsiN(3.xXM(6)) 

P  =  XM(l)x(l .-XM(2)xx2. ) 

RR  =  P/(1.+XM(2)XC0S(F)) 

XSPA=  J(2)XRXR/XM(11X((XM<1)/  RR  )XX3 . X( ( 1 . -3 ./2 . XSINCXMC 3) )xSIN 
1(XM(3)))+3./2.xSIN(XM(3))xSINCXM(3))xC0S(2.x(XMC4)+F)))- 
2(l.-3./'2.XSIN(XM(3nxSIN(XM<3n)x(l.-XM(2)xx2.)xx(-l.5D0)) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
XX  XX 

XX  FUNCTION  XSPE  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


FUNCTION  XSPE(XM) 

DIMENSION  XM(6) 

REAL  J ( 4 ) 

DATA  J,R/1 . ,  0 .001 08276  DO, -0 . 0  00002S5D0, -0 . 00000156  DO, 6  378 . 165D0/ 
r  =  XH(6)+(2.XXMC2)-XMC2)XX3./4+5./96.XXM(2)XX5. )XSIN(XM(6)) 

1+(1 .25D0XXM(2)XX2.-11 ,/24 . XXMC 2 ) xx4 . + 17 ./192 . XXMC2 >xx6 . ) xSIN( 2 . x 
2XMC 

36))+<13./12.xXM(2)xx2.-43./192.xXM(2)xX5. )XSIN(3.XXM(6)) 

P  =  XM(1)X(1 ,-XM(2)xx2. ) 

XSPE  =  J(2)/2.x(R/P)xx2.x<l.-1.5D0XSIN(XM(3))xSIN(XM(3))) 
1X((1.+1.5DOXXM(2)XXMC2)-(1.-XM(2)XXM(2))XX1,5DO)/XM(2)+3.X(1 .+ 
2XMC21 

3XXfU2)/4. )XC0S(F)  +  1 ,SDOXXM(2)XCOS(2.XF)+XM(2)XXM(2)/4.XCOS(3.XF) ) 
4+3./8.xj(2)x(R/P)xx2.xSINCXM(3))xSIN(XMC3))x((l ,+H ,/4.xxM(2)xx2. ) 

5  *C0S(2.xXM(4)+F)+XMC2)xXM(2)/4.xC0S(2.xXM(4)-F)+5 .*XM(2)x 

6  C05<2.XXM(4)+FX2. )+(7 .+17 ,/4.*XM(2)xxM(2))/3.xC0S<2.XXM(4)+3.*F) 

7  +  1 . 5D0XXMC  21 *C0S(2 . xXM(4 )+4 . xF)+XM( 2)XXM( 2 )/4 . XCOSL  2 . XXM<  4)  +  5 . xF) 
3+1 ,5DOXXM(2)XCOS(2.XXM(4)) ) 
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RETURN 

END 

x 

X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

x  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  xx  FUNCTION  XSPI  XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

FUNCTION  XSPI(XM) 

DIMENSION  XMC6 ) 

REAL  J ( 4 ) 

DATA  J>R/1. >0.001 0827 6 DO >-0.0000025500,-0. 000001 56  DO >6378.16500/ 
F  =  XMC6)+(2.*XM(2)-XMC2)*x3./4+5./96.*XM(2)**5.)*SIN(XM<6)) 

1+C1 .25D0XXMC2) **2.-11 ./24 . XXMC 2)**4 . +17 . /I 92 . *XM( 2) *X6 . )XSIN(2.* 
2XM( 

36 ))+C13./12.*XM< 2) **2. -43./192. *XMC2)**5. )*SIN( 3. *XM(6)) 

P  =  XM(l)x<l ,-XM(2)xx2. ) 

XSPI  =  3./8 .XJ(2)*(R/P)**2.XSIN(2.XXMC3))*CXM(2)*C0S(2.XXMC4)+F) 
1  +C0S(2.X(XM(4)+F))+XM(2V3.XC0S(2.XXM(4)  +  3.XF)) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

XX  FUNCTION  XSPW  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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X 
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X 

X 

X 

X 
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FUNCTION  XSPWCXMl 
DIMENSION  XM(6 ) 

REAL  J<4) 

DATA  J>R/1., 0.0010827 6  DO >-0.00000255 DO >-0.000001 56 DO >6 37 8. 165D0/ 

F  *  XM(6)+C2.XXM(2)-XM<2)XX3./4+5./96.XXM(2)*x5. )xSINCXM(6>) 
1+(1.25D0XXM(2)XX2.-11 ./24 . XXMC2)**4 . +17 ./192 . XXM(2)XX6 . )*SINC2 . * 
2XM( 

26))+<13./12.xXM(2)xx2.-43./192.*XM(2)xx5. )XSIN(3.XXM(6)) 

P  =  XMCl)x(l .-XM(2)**2. ) 

XSPW  =  3./4.xJC2)*<R/P)**2.x(4.-5.*SIN(XMC3))*SIN<XMC3))) 
1*(F-XM(6)+XM(2)*SIN<F))+1 .5D0XJC2)*(R/P>**2.*(1.-1 .5D0xSIN(XMC3)) 
2**2.  ) 

3*((l.-XMC2)*XM(2)/4. )/XM< 2) *SIMC F)+ . 5D0*SINC 2 . *F)+XMC2) XSINC 3 . *F) 
4/12 

5.  )-l .5D0XJ(2)*(R/P)**2.*<(SINCXM(3) )*SIN(XM(3))/4.+XM(2)**2./2.*( 

61  .-15./8.*SIN(XM(3))**2. ) )/XM( 2) *SIN( 2 . *XM( 4 )+F)+XM( 2)/ 16 .*SIN 
7(XM(3))**2.XSIN(2.*XM(4)-F)+.5D0*(1 . -2 . 5D0*SIN( XMC 3 ) ) **Z . )*$IN(2.* 

8(XM(4)+F))-(7./12.*SIN(XM(3))**2.-XMC2)*XI1(2)/6.*(l.-19./8.*3IN 
9CXMC3) )**2. ) )/XM( 2)*SIN( 2 . *XM( 4 )+3 . *F)-3 ./3 . *SI N( XM( 3 ) ) **2 . *SIN 
11(2.XXM(4)+4.XF)-XMC2)/16 .XSIH(XM(3) )**2.*SIN(2.*XM(4)+5.*F)) 
12-9./16.*J(2)X(R/P)X*2.*SIN(XM(3))**2.*SIN(2. *XM<  4) ) 

RETURN 

END 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

**  FUNCTION  XSPO  *x 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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FUNCTION  XSPO(XM) 

DIMENSION  XMC6 ) 

REAL  JC4) 

DATA  J.R/l . ,0.00108276  00,-0 . 00000255D0, -0 . 00000156 DO, 6 378 .16  5D0/ 

F  =  XMC6)+C2.*XMC2)-XMC2)**3./4+5./96 ,*XMC2)**5. )*SINCXMC6)> 

1+C1 .25D0*XMC2)**2.-11 ./24.*XMC2)**4.+17 ./I 92 . XXM( 2 )XX6 . )xSIN(2.x 
2XM( 

36))+C13./12.*XMC2)**2.-43./192.*XMC2)**5. ) XSINC 3 . XXM( 6 ) ) 

P  =  XMC 1 )*C 1 . -XM( 2 )*x2 . ) 

XSP0=-1 . 5D0*JC2)*CR/P)**2.*C0SCXMC3))*CF-XMC6)+XMC2)*SINCF)-XMC2) 
l/2.*SINC2.*XMC4)+F)-SINC2.*CXMC4)+F))/2.-XMC2)/6 .*SINC2.*XMC4)+3.x 
2F) ) 

RETURN 

END 


KXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
XX  XX 
xx  FUNCTION  XSPM  xx 
XX  XX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


FUNCTION  XSPMCXM) 

DIMENSION  XM(6 ) 

REAL  J(4) 

DATA  J.R/l. ,0.001 08276  DO, -0.0000025500,-0. 0000  015600,6378.16500/ 

F  =  XMC6)+C2.*XMC21-XMC2)**3./4+5./96.*XMC2)**5.)*SINCXMC6)) 
1+C1.25D0*XMC2)**2.-11 ./24 . XXM( 2)**4 . +17 ./ 192 . *XMC2)**6 . )*SINC2.» 
2XMC 

36mC13./12.*XMC2)**2.-43./192.*XMC2)X*5.  )*SINC 3 . *XMC6 ) ) 

P  =  XMC 1 ) *C 1 .  -XMC  2) **2  .  ) 

XSPM  =  -1.5D0XJ(21X(R/P)XX2.XSQRT(1.-XM(27XX2.)/XM(2)X((1.-1.5D0X 
1SIN( 

2XMC3»xx2.)x(Cl.-.25D0*XM<2)*XMC2))*SINCF)+XM(2)/2.*SINC2.*F)+ 

3XMC2) 

4*XMC2)/12.*SIN<3.*F))+.5D0*SINCXM(3) )*SINCXMC3) )*CC-.5D0)*C1 .♦ 
51.2500 

6*XMC21XXM(2n*SINC2.*XM(4)+F)-XM(2)*XMC2)/8.*SINC2.*XMC4)-F)+7  ./6  . 
7*Cl.-XMC2)**2./28 . )*SINC2 . *XMC 4 )+3 . *F)  + . 75D0XXMC 2) »SIN( 2 . XXM( 41+4 . 
8*F 

9)+XM(2)xx2./8 .xSIN(2.xXM(4)+5.xF)))+9./16 .XJ(2)X(R/P)xx2.x 
1SQRTC1.-XMC2)**2. )*SIN(XMC3) ) XSINC XMC 3 ) ) XSINC 2 . *XMC 4) ) 

RETURN 

END 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
XX  XX 
XX  FUNCTION  XSPR  ** 
XX  XX 


X 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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X 

X 

X 
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FUNCTION  XSPRCXM) 

DIMENSION  XMC  6 ) 

REAL  JC4J 

DATA  J,R/1  . ,  0.0010827600,-0 . 00000255DO , -0 . 000001 56 DO , 6378 . 16  5D0/ 

F  =  XMC6)+C2.*XMC2)-XMC2)**3./4+5./96 .*XMC2)**5. )*SINCXMC6)) 

1  +  C1 .25D0XXM(2)XX2.-11 . /24 . XXMC 2 ) **4 . +17 . / 1 92 . *XMC 2) **6 . ) XSINC 2.x 
2XMC61 )+(13./12.xXM(2) **2.-43. / 192. *XMC2)xx5. ) XSINC 3. XXMC6) ) 

P  =  XMC11XC1 .-XMC2)**2.  ) 

RR  =  P/Cl  .+XMC2)*C0SCF) ) 

XSPR  =  -.5D0*JC2)*CR**2./P)*C1.-1.5D0*SINCXMC3))*SINCXMC3)))*C1.+ 
1XMC2) 

2/ Cl ,+SQRTCl .-XMC 2)**2. ) ) *C0SC F)+2 . / C SQRTC 1 . -XMC2)**2. ))*RR/XMC1 
3))+.25*JC21*R**2./P*SINCXMC3))**2.*C0SC2.*XMC4)+2.*F) 


*  *  *  *  K 


RETURN 

END 

x 

X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

*  XX  SUBROUTINE  MEAN  xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

SUBROUTINE  MEANCXOSC,XM) 

COMPUTES  THE  MEAN  ORBITAL  ELEMENTS  FROM  THE  OSCILATORY  ELEMENTS 
AND  THE  FIRST  ORDER  PERIODIC  VARIATIONS  BY  THE  RELATIONS, 

XM  *  XOSC  -XSPCXM) . 


DIMENSION  X0SC(6),XM(6),XMN(6 ) , XSP( 6 ) , TOL ( 6 ) , DELXM( 6 ) 

DATA  TOL/. 01  DO, . 00001  DO,. 00001  DO,. OOOOIDO,. 00001  DO,. 00001  DO/ 

KOUNT  =  0 

X...  INITIAL  GUESS  FOR  MEAN  ELEMENTA. 

DO  10  1=1,6 
10  XMC I )  =  XOSC(I) 

X...  CALCULATE  THE  SHORT  PERIODIC  VARIATIONS. 

20  XSPC1)  =  XSPA(XM) 

XSP(Z)  =  XSPECXM) 

XSP(3)  =  XSPI(XM) 

XSPC4)  =  XSPW(XM) 

XSP(5)  =  XSPO(XM) 

XSP(6 )  =  XSPM(XM) 

KOUNT  =  KOUNT+1 
DO  50  J=l,6 

XMN(J)  =  XOSCCJ)-XSP(J) 

DELXMC  J )  =  ABS( XMN( J )-XM( J ) ) 

IFCDELXM(J) .GT.TOLCJ))  XM(J)  »  XMN(J) 

50  CONTINUE 
DO  60  K=1 , 6 

IFCDELXM(K).GT. TOLU). AND. KOUNT. LT. 30)  GO  TO  20 
60  CONTINUE 

70  IFC KOUNT . GT . 30 )  WRITE<6,100)  DELXM 
90  CONTINUE 

100  F0RMATC/T5, 'CONVERGENCE  OF  ONE  OR  MORE  MEAN  ELEMENTS  WAS  NOT  MET.' 
1/ 'THE  ABSOLUTE  DIFFERENCE  IN  THE  ELEMENTS  ARE  SHOWN . '//7E12 . 5 ) 
RETURN 
END 
X 
x 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX  SUBROUTINE  GEOP  xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  GEOP  (XM.DXDT) 

DIMENSION  XM( 6 )  ,  DXDT ( 6 ) 

REAL  N , MCONS , NCONST , NCON 
REAL  J ( 4 ) 

DATA  J,R/1. ,0.00!08276DO,-0.00000255DO,-0.00000156DO,6378. 165D0/ 
DATA  U/398600 ./ 

N  =  SQRT(U/(XM(1)xx3. ) ) 

P  =  XM<1)X(1 ,-XM<2)xx2.  ) 

DXDT ( 1 )  =  0. 

S3  =  SIN(XM(3))XSIN(XM<3)) 

XS  =  1.-XMC2)XX2. 

CS  =  COS<  XMC  3 ) ) 

DED  =  (-3. )/32.XNxj(2)xj(2)X(R/P)xx4.xS3XC14.-15.xS3)xxM(2)x 

1XSXSIN(2.XXM(4)) 
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DEDO  =  3./3.XNXJ(3)x(R/P)xx3 . XSINC XMC 3 ) ) x( 4 . -5 . xS3 )XXSXC0S( XMC 4 ) ) 
DEDOT  =  15./32.XNXJ(4)X(R/P)XX4.XS3X(6.-7.XS3)X(XM(2)XXS) 
lxSINC 2 . XXMC 4 ) ) 

DXDT ( 2 )  =  DED-DEDO-DEDOT 

DID  =  3./64.xNxj(2)xj(2)x(R/P)xx4.XSINCXMC3)x2. )x(14.-15.xS3)x 
1XM(2)XX2.XSIN(2. 3XXM  C  4 ) ) 

DIDO  =  3./8.xfl*J(3)x(R/P)**3 ■ XCOSC XM( 3 ) )x( 4 . -5 . XS3)XXM( 2) XCOSC XMC 4 
1) ) 

DIDOT  =  15./64.xNXJ(4)x(R/P)xx4.xSIN(2.xXMC3))x(6.-7 ,XS3)XXMC2)X 
1XM(2JXSIN(2.XXM(4) ) 

DXDT  C  3 )  =  DID+DIDO+DIDOT 

DWD  =  .75D0XNXJ(2)X(R/P)*X2.X(4.-5.XS3) 

DU  =  2.X(14.-15.XS3)XS3-C2S.-153.XS3+I35.XS3XX2. )XXM(2)XX2. 

DWDO  =  3 ./16 . xnxJ(2)x*2 . X(R/P)xx4. x( 48-103 . XS3  +  215 . /4 .  XS3XX2  .  + 
1(7  .-4.5D0XS3-45./S .XS3XX2. )xxM(2)xx2.+6 ,x(l ,-1.5DQxS3)x(4.-5.xS3)x 
2  SQRT(XS)-. 25DQXDWXC0S ( 2 . XXMC  4 ) ) ) 

X...  PRINT  X,  XMC  2) ,  XMC  3 ) ,  XMC 4 ) . R, P , N, J ( 3 ) , S3 

DWDOT  =  3./8 .XNXJ(3)X(R/P)XX3X( C4 . -5 . XS3 ) x( S3-XMC 2 )XX2XC0S CXMC 3 ) 

1) XX2)/CXM(2)XSIN(XM(3)))+2.XSINCXMC3))X(13.-15.XS3)XXMC2))XSIN(XM 
2(4) ) 

DUN  =  S3»(6 . -7 . XS3)-C6 . -35 . X53+31 . 5x$3X*2 . )XXM( 2)X*2 . 

D1'ID0TT  =  15./32.XNXJ(4)X(R/-P)X^4.XC16.-62.XS3+49.XS3XX2.  +  .75D0X(24.- 
184.XS3+63.XS3XX2. )XXM(2)XX2.+DWWXC0SC2.XXMC4))) 

DXDTC4)  =  DWD+DWOO+DWDCT-DWDQTT 
DOD  =(-l .5D0)xnxj(2)X(R/P)xx2. XCS 

DODO  =  1 .5D0XNxj(2)xx2.x(R/P)xx4.xCSXC2.25D0+1 .5XSQRTCXS)-S3 
1XC2.5D0+2.25D0XSQRTC  XS ) ) +XMC2)XX2 . /4 . XC 1 . +1 . 25D0XS3)-XM( 2) XX2 . 
2/8.x(7.-15.xS3)xC0SC2.xxMC4))) 

DODOT  =  3./8 .XNXJ(3)X(R/P)*X3.X(15.XS3-4)XXM(2)XCS/SIN(XM(3) )xSIN( 
lXMC  4)  ) 

D0D0TT=15./16 .XNXJ(4)X(R/P)XX4.XCSX( C4.-7.XS3)XC1 .+1 .5D0XXMC2)XX2. 

2)  -(3. -7 .X$3)XXM(2)X*2.XC0S(XM(4)X2. )) 

DXDT ( 5 )  =  DOD-DODO-DODOT+DODOTT 

NCON  =  3./8 .XNXJ(2)X*2. 

NCONS  =  NXJ(4)X(R/P)XX4. 

NCONST  =  3./8 .XNXJ(3)X(R/P)X*3. 

DM  =  Nxa.+l.5D0XJ(2)x(R/P)xx2.XCl.-1.5D0xS3)xSQRTCXS)) 

DMD=1 .5D0xNxJ(2)X*2x(R/P)xX4X(C1 .-1 ,5D0xS3)xx2xXS+(l .25D0XC1 .- 
12.5D0X 

1S3+13./8.xS3XS3)+5./8.X(1.-S3+5./8.xS3XS3)XXMC2)x*2+S3/16.x14.- 
2  15.XS3)X(1 .-2.5D0XXMC2)XX2)XC0SC2.XXMC4)))XSQRT(XS)) 

DM1 =<3.-7 .5D0XS3+47 ,/8 . XS3XS3+C 1 . 5D0-5 . XS3+1 17 . /6 . XS3XS3 ) XXMC 2 ) x*2 
1- C 1 .+5.XS3-101./8.XS3XS3)/8.XXMC2)XX4.  ) 

DM2  =  XMC2)XX2./8.XS3X(70.-123.XS3+C56 .-66.XS3)XXMC2)XX2. )XC0SC2.X 
lXMC  4 ) ) +27 . / 1 28 . XXM(2)XX4 . XS 3XS3*COS  C  4 . XXMC  4 ) ) 

DMDO  =  NC0Nx(S/P)xx4 ./(SQRTCXS) )x( 3 . XDM1+DM2) 

DMDOT  =  NCONST XS IN (XMC  3) ) *<  4 . -5.*53)XC1 . -4. XXMC  2)XX2. )/XM(2)x 
1SQRT ( XS ) * 

lSIN(XM(4))-45./128.XNC0NSxC8.-40.xS3+35.XS3xS3)xXMC2)xx2.xSQRT(XS) 
DMD0TT=15./64.XNC0NSXS3X(6 .-7.XS3)X(2. -5.XXMC2)XX2. )xSQRT(XS)xCOSC 
12 . XXMC  4 ) ) 

DXD~(6)  =  DM+DMD+DMDC-DMDOT+DMDOTT 
RETURN 
tMD 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XX  XX 

XX  SUBROUTINE  SOLOC  XX 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  SOLOC  (  TVE  ,  T  ,  BLA  ,  DS  ,  AS  ,  AB  ) 

DATA  EPS  1/  23.5D0  /  ,  PI  /  3.14159D0/,  TOL/  0.0001D0  / 
PI1  =  0.5DOXPI 
P 1 2  =  2.  x  PI 
PI 3  =  1 . 5 D 0 x  PI 


x 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 
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XHKXlXlXXHtX  5 *********** 


RPD  *  PI  /  180. ODO 

DELT  »  T  -  TVE 

EPS  *  EPS1  *  RPD 

CLON  =  0 . 985647 DO  *  RPD  *  DELT 

IF  (CLON . LT .  0.0  )  CLON  =  PI2  +  CLON 

IFLCL0N.GT.PI2)  CLOU  *  AMODCCLON, PI2) 

SDS  *  SIN(EPS)  x  SIN  (  CLON  ) 

DS  *  ASIN  (  SDS  ) 


8LAR  »  BLA  x  RPD 
>0  10  I  ®  1  ,  5 
J  =  I  -  1 

ARG  a  0.5D0  *  J  *  PI 
TST  =  ARG  -  CLON 
TST1  »  AB$<  TST  ) 

IF  (  TST1  .  LE  .  TOL  )  GO  TO  20 
10  CONTINUE 

ETA1  =  TAN(DS)  /  TAN(EPS) 

ETA  =  ASIN(ETAl) 

ETA  a  ABS(ETA) 

IF  ((  CLON  .  LT  .  PI1  )  . 


IF  ((  CLON  .  LT  .  PI1  )  .  AND  .  (  CLON  .  GT  .  0.0))  AS  *  ETA 
IFCCCLON.LT. PI) .AND. (CLON. GT. PI1))  AS  =  PI  -  ETA 
IF  ((CLON. LT.PI3). AND. (CLON. GT. PI))  AS  =  PI  +  ETA 
IF  CCCL0N.LT.PI2). AND. (CLON. GT.PI3))  AS  *  PI2  -  ETA 


GO  TO  21 
20  AS  *  CLON 


21  A3  =  AS  +  BLAR 
RETURN 
END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XXXKXXXXXXXXXXXXXXXXXXXXXXXXXX 
XX  XX 

XX  SUBROUTINE  DATPRP  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


SUBROUTINE  DATPRP  (  X,Y,XX,YY  ) 

COMMON  /  PREP  /  ASUN,  DSUN,  RAS,  DECS 

X  »  ASUN 

Y  a  DSUN 

XX  «  RAS 

YY  =  DECS 

RETURN 

END 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

xx  SUBROUTINE  PERIOD  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  PERIOD  (  AO,  ES,  XI,  PA  ) 

REAL  J2 

COMMON  /  XXXX  /  IET 

DATA  J2  /  0.00108276D0/,  R/  6378.165D0/  ,  U/3986 00 . / . PI/3 . 1 A 159D0/ 
PA  a  2.XPIX(SQRT(A0XA0XA0/U) )X(1 .-(3.XJ2XRXRX(3.XSIN(XI)XSIN(XI) 

1  -2. 0)/(4.XAOXAOX((1 .-ESXESJXXI ,5D0)) ) ) 

IF  (  IET  .  E<3.  0  )  RETURN 
WRITE  (6,100)  PA 

100  FORMATC/lOX, 'INITIAL  ANOMALISTIC  PERIOD®  ',F12.3,'  SEC*,//) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


(A 


*  *  *  K  *  *  *  tt  *  *  *  *  )*  *  J*  *  !K  *  K 


R 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
XX  XX 
XX  SUBROUTINE  SALT  xx 
xx  xx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  SALT  ( D, XM, RM, Z, RE  ) 

DIMENSION  XMC6) 

DATA  ROE/6373 . 165D0/, E0E2/6 . 69342E-03/ 

RE  =  ROE/S QRTCl.+C  E0E2/C l.-E0E2)J*SIN(D)xSIN(D)) 

Z  =  RM  +  XSPRCXM)  -  RE 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
XX  XX 
xx  SUBROUTINE  SATLOC  xx 
XX  XX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  SATLOC  CMTYPE, XM, DELT, DSUN, ASUN, DECS, RAS, GMGLT , R) 
DIMENSION  XM(6 ) 

DATA  B0/398603. 254DQ/,W/.7292115£-04/, PI/3. 14159D0/ 

AM  =  0.0 
PI02  =  0.5D0XPI 
PI2  *  2 . *PI 
P2  *  PI2 
PI03  =  3.XPI/2. 

AG  »  AMOD  (  WXDELT  ,  PI2  ) 

A  =  XMC 1 ) 

CE  *  XM( 2) 

Cl  *  XM( 3) 

G  =  AMOD  C  XMC 4 )  ,  PI2  ) 

H  »  AMOD  (  XMC  5 )  ,  PI2  ) 

CALL  POSVELCA,CE.CI,AM,G,H,BO,R1,R2,R3,V1,V2,V3,R,V) 

DECS  *  ATAN  C  R3/SQRT C  RlXRI  +  R2*R2  )) 

RAS  *  ATAN2  (  R2  ,  R1  ) 

IF  C  MTYPE  .EQ.  0  )  GO  TO  10 

ALP  *  TANCDSUN)  x  COSCCI)  /  SINCCI) 


ALP  *  ASINCALP) 
ALP  =  ABSCALP) 
DECS  =  DSUN 
IF  C  MTYPE  .  EO 


IF  C  MTYPE  .  EO  .  1  )  DECS  =  -  DSUN 
U1  =  ASINC  SINCDSUN)  /  SINCXMC3) ) ) 

U2  =  ASIN  C  SIN(DEC5)/SINCXMC3) ) ) 

T1  *  ASUN  +  PI02 

T2  »  ASUN  -  PI02 

HU  a  g 

OMA  =  H 

IF  CT2  .LT.  0.0  )  GO  TO  20 
IP  C  T1  .GT.  PI2  )  GO  TO  30 

IF  (COMA  ,LE.  Tl)  .AND.  COMA  .GE.  T2))  GO  TO  40 
GO  TO  50 

20  T 3  =  PI2  +  T2 

IF  C ( ( OMA . L  E . Tl ) . AND . C  OMA . GE . 0 . 0 ) > . OR . C  C  OMA . L  E . P2 ) . AND . C  OMA . GE . T 3 
1  )))  GO  TO  40 

GO  TO  50 

30  T3  »  Tl  -  PI2 

IFCCC0MA.GE.T2) . AND . C OMA . LE . PI  2 ) ) . OR . C C OMA . GE . 0 . ) . AND . C OMA . L E . T3 
1  )))  GO  TO  40 

GO  TO  50 

40  IF  C  Cl  .  GE  .  PI02  )  GO  TO  41 
RAS  =  OMA  +  ALP 

IF  C  DECS  .  LT  .  0.0  )  RAS  *  OMA  -  ALP 
GO  TO  42 

41  RAS  =  QMA  -  ALP 


B 


% 


rr)T>:nTii 


*  *  IK  *  *  *  *  5X  K*  *  *  *  *  *  *  *  *  *  *  * 


nranos 


aiim  u  MJawnuminurxy  xjrtrxy  Ryry^jrorarx  ray*  ywyar".  y»r*y*>r».>rvy  *x  -.x-j,  -r  ->  -. 


IF  (  DECS  .  LT  .  0.0  )  RAS  =  OMA  +  ALP 
42  IF  C  MTYRE  .  EQ  .  1  )  RAS  *  RAS  +  PI 
RAS  =  AMOD  (  RAS  ,  PI2  ) 

GO  TO  10 

50  IF  (  Cl  .  GE  .  PI02  )  GO  TO  51 
RAS  «  OMA  +  PI  -  ALP 

IF  (DECS  .  LT  .  0.0  )  RAS  =  OMA  +  PI  +  ALP 
GO  TO  42 

51  RAS  =  OMA  +  PI  +  ALP 

IF  (  DECS  .  LT  .  0.0  )  RAS  *  OMA  +  PI  -  ALP 
GO  TO  42 

10  ELNG  *  RAS  -  AG 

IF  (ELNG  .LT.  0.0  )  ELNG  =  PI2  +  ELNG 
ELNG  =  AMOD(ELNG, PI2) 

GMGLT  =  ASIN  (  0 . 9792D0*SIN(DECS)+0 . 2028D0XCOS( DECS) *COS( ELNG  - 
1  ( 291 . XPI/180 .  )  ) ) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

XX  SUBROUTINE  DRAG  XX 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  DRAG(  XM,  DEL.  TT,  DELMG  ) 

REAL  N. 10, II, 12, 13, 14, 15, 16, 17 

COMMON  /  SOLAR  /  F107, FBAR, AKP.TVE 

COMMON  /  PREP  /  ASUN.DSUN, RAS, DECS 

DIMENSION  XM(6 ) ,  SUN(4)  ,  SAT(2)  ,  GE0(3>  ,  BF(8) 

DATA  B0/398603.254D0/,BLA/30./,PI/3. 14159D0/ , EOE/ . 00335D0/ 

DATA  IYR/1981D0/,TVD/1./,THR/17./,TM/14./,TS/56.243D0/ 

PI2  =  2.  x  PI 

SUN(3)  =  23. 5D0X(p 1/180.) 

SUN( 4)  *  1.0 

GEO(l)  »  F107 

GE0(2)  *  FBAR 

GEO( 3 )  *  AKP 

N  =  SQRT(B0/(XM(1)XX3)) 

TCOR  »  AMOD  (  XM( 6 )  ,  PI2  ) 

TCOR  »  TCOR  /  N 

TTT  »  TT  -  (  TCOR  /  86400.  ) 

TTR  *  1721044.  +  367XIYR  -  (7XIYR/4)  +  TVD  +  ((  THRX3600.  +  TMX60. 
I  +  TS  )  /  86400.  )  -  2400001. ODO 
DELT  *  TTT  -  TTR 

CALL  SOLOC(TVE, TTT, BLA, DSUN, ASUN , ABLG ) 

MTYPE  =  0 

CALL  $ATLOC(MTYPE, XM, DELT, DSUN, ASUN, DECS, RAS, GMGLT, RM) 

3TR  =  XM(6 ) 

XM(  6  )  a  0.00 

CALL  SALT(  DECS, XM, RM, HP, RE  ) 

XM( 6 )  =  STR 

CALL  DATPRP(  SUN( 1 ) , SUN(2) , SAT( 1 ) , SAT( 2)  ) 

CALL  ATMDEN( TTT, SUN, SAT, GMGLT, HP, GEO, RHO,DRDH,SH) 

SHOP  =  RHO  X  (  10.  XX9.  ) 

S HP  =  SH  X  0 . 001  DO 

CALL  JAC60( DECS, DSUN, RAS , ASUN, BL A , HP, FI 07 , RHOP , SHP) 

MTYPE  '  1 

CALL  SATLOCCMTYPE, XM, DELT, DSUN, ASUN, DECS, RAS , GMGLT, RM) 

CALL  DATPRP  (  SUN( 1 ) , SUN( 2 ) , SAT( 1 ) , SAT ( 2 )  ) 

CALL  ATMDEN(TTT, SUN, SAT, GMGLT, HP, GEO, RHO, DRDH.SH) 

RHOMIN  a  rho  x  (  10.XX9  ) 

SHMAX  a  sh  X  0.001D0 

CALL  JAC6  0( DECS, DSUN, RAS, ASUN, BLA, HP, FI 07, RHOMIN, SHMAX) 

MTYPE  =  2 

CALL  S ATLOC( MTYPE, XM, DELT , DSUN, ASUN, DECS , RAS , GMGLT, RM) 

CALL  DATPRP  (  SUN(l),  SUNC2),  SAT(l),  SATC2)  ) 

CALL  ATMDEN(TTT, SUN, SAT, GMGLT, HP, GEO, RHO, DRDH,SH  ) 
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$ 


a 


X  RHOMAX=RHO  X  (10.XX9  ) 

X  SHMIN  =  SHX0.001D0 

CAL  .  JAC60 (  DECS , DSUN, RAS, ASUN, BLA, HP, F107 , RHOMAX, SHMIN) 

A  *  $IN(DSUN)XSIN(XM(3))xSIN(XM(4) )+C3S( DECS )x( COS( XM( 5)-ABLG)x 
1  C0S(XM(4))-C3S(XM(3))XSIN(XM(5)-ABLG)XSIN(XM(4)) ) 

B  =  SIN(DSUN)x$IN(XM(3))xC0S(XM(4))-C0S( DECS )X( COS( XM( 5)-ABLG) x 
1  SIN(XM(4))+C0S(XM(3) ) XSIN( XM( 5 ) -ABLG) XC0S( XM( 4) ) ) 

F  =  (RHOMAX  -  RHOMIN)  /  (RHOMAX  +  RHOMIN  ) 

RHOO  =  RHOP  -  0.5D0X(RHOMAX  +  RHOMIN)  x  F  X  A 

Z  *  XMC 1 )  X  XM( 2)  /  SHP 

Z1  =  ZX(2  -  1.)  -  (  FXA  /  (1.  +  F  X  A  )) 

Z2  =  (  1.  /  (  1.  +  F  X  A  ))  -  Z  X  (  1.  -  0.5D0XZ)  -  0.5D0 
TM1  =  -Z1XZ1  -  2.  X  ZXZXZ2 
CTHAVG  =0.0 

IF  (  TM1  .GT.  0.00  )  CTHAVG*  (Z1  +  SQRKTM1 ) )/( ZXZ) 

SHNM  =  ABS(  SHMAX  -  SHMIN  ) 

HAVG  =  0.5D0X(SHMAX+SHMIN)+  (  SHNM  /  (SHMAX+SHMIN) )XAXCTHAVG 
BETA  =  1.0  ✓  HAVG 

C  *  0.5DOXBETAXEO£XSIN(XM(3))XSIN(XM(3))X(RE+HP) 

CC  *  C  x  C 

ZZ  *  BETA  X  XM(1 )  X  XM( 2) 

ALP  *  0.00 
MO  =  1 
NN  =  8 

CALL  MMBSIR(ZZ,ALP,NN,MO,BF.NZ) 

10  =  BF( 1 ) 

11  =  BF(2) 

12  =  BF( 3 ) 

13  *  BF( 4) 

14  *  BF(5) 

15  *  BF(6 ) 

16  *  BF(7 ) 

17  *  BF(8) 

FA  *  FXA 
FB  *  FXB 

E  *  XM( 2) 

C2W  =  COS( 2 . XXM( 4) ) 

S2M  *  SIN( 2 . XXM( 4) ) 

C4W  =  COS( 4 . XXM( 4) ) 

S4W  =  SIN(4 . XXM(4) ) 

DELMG=-DELX(SQRT(BOXXM(I)))XRHOOXEXP(-ZZ-CXC2W)X((1.+0.25DOX 

1CXC)X 

1(I0+EXI1)+FAX((1.+0.25D0XCC)X(I1+EXI2))+0.5D0XCX((2.XI2-EX(I1-3.X 
213)  )  +  FAx( H  +  I3+2. XEXI4) )xC2W-0 .5D0XCXFBX(( H-I3)+2.XEX( 12-14) )x 

3  52W  +0.125D0XCXCX((2.XI4-EX(3.XI3-5.XI5) )+FAx( ( 13+15 )+Ex( 3 . XI6-I2 

4  )))XC4W  +  0.125D0XCXCXFBX((I5-I3)+EX(I2-4.XI4+3.XI6))XS4W) 

XT1  =  DEL x( SORT ( B0XXM( 1 ) ) )XRHOO 

XT2  *  DELMG  /  XTI 

XT3  *  EXP(-ZZ)X(I0  +  EXI1) 

XT4  =  XT 2  /  XT3 
DELMG  =  1000.  x  DELMG  /  9.8D0 
X...  PRINT  x,  ZZ,  RHOO,  XTI,  XT2,  XT 3,  XT 4,  DELMG  ,  DEL 
RETURN 
END 
x 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX  SUBROUTINE  JAC60  xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

SUBROUTINE  JAC60  ( D, DS , A , AS , BLA, Z, FI 07 , RHO, SH  ) 

BL  =  BLAX3 . 14159D0  /  180. 

XP  *  -16.021D0  -  O.OQ1985DOXZ  +  6.363D0xEXP(  -0.0026D0XZ) 

RO  *  10.XXXP 

CSPSI  =  SIN( D)*SIN( DS)  +  COS( D)xCOS( DS ) XC0S( A-AS-BL ) 

F  =  (0.5D0XC1.  +  CSPSI ) )XX3  . 

RHOl  *  1.  +  0 . 19D0X( EXP( 0 .0055D0XZ)  -  1.9D0)  XF 
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M)K)K3K)K)K)K2K)K)K)K 


RHO  *  O.OIDQ  X  F107  <  RO  (  RHOl  X  (10.XX12) 

DRODZ  »  2.3026DO*ROX(-.00198SDO  -  . 016544D0XEXPC - . 0026D0XZ) ) 
DR1DZ  =  .001045D0XFXEXPC .0055D0XZ) 

SH  =  -CROXRHOIVCROXDRIDZ+RHOIXDRODZ) 

RETURN 

END 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
XX  *  XX 

xx  SUBROUTINE  THRST  xx 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  THRST (  XM,  DXM,  IRC,  K,  PA,  DWOA  ) 

DIMENSION  XM( 6 )  ,  DXM(6),  IRV(IO),  DAC1Q) 

COMMON  /  OADJ  /  NOA,  IRV,  DA,  SPI,  WT,  PSENS 
DATA  BO  /  398603. 254D0/,  PI/3 . 14159D0/ 

DO  10  I  »  1,6 
DXM(I)  =0.0 
10  CONTINUE 
DWOA  =  0.0 

IF  (  IRC  .EQ.  0)  RETURN 

THT  =  XM(6)  +  2.XXMC2)XSIN(XMC6))  +  (5 ./4 . )XXM(2)XXM(2)XSINC2 . XXM 
1  (6)) 

R  »  XM(1)X(1 .  -  XM(2)XXM(2»/U.  +  XM(2)XC0S(THT) ) 

V2  =  BO  x  (C2./R)  -  <1./XM(1))) 

V  *  SQRTCV2) 

DVI  =  DA(K)xB0/(2.XXMaiXXM<l)XV) 

DXM( 1 )=  DA(K) 

DXM(2)  =  2.X((C0S<THT)+XMC2))/V)XDVI 
DXMC4)  =  2.XSIN(THT)XDVI/CXM(2)XV) 

DXM<6 )  =  (2.XPI/PA)-U./(XM(1)X(SQRTC1.-XMC2)XXMC2)))XV))X2.X 
1  SIN(THT)X<  (XM(  1  )x(  1 . -XM( 2) *XM( 2)  l/XMC 21 )  +  RXXMC2DXDVI 
PA  =  3 . XPAXXMC 1)X(V/B0)XDVI 
DWOA  =  (WT/(SPIX9.8D01)X1000.*DVI 
IF  (PSENS  .EQ.  0)  THEN 

PRINT  x,  THT ,R,V,IRV(K),DA(K),DVI, DXMC2) , DXM( 4) , DXM(6 ) , PA, 

1  DWOA 

END  IF 
RETURN 


!**m*m*o 
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APPENDIX  B 

SAMPLE  INPUT  AND  OUTPUT 


05000, 
2.2,  . 
100., 

1  10947 

2  10947 

3  10947 

4  10947 

5  10947 
CCCC  X. 
CD  1. 
A  0. 


001,0,1,0,025 

0000500,  1.0,  230. ,100. ,09100. 
100.,  2.0  ,43222.7382 
U  90070  78  060  A  0 


1.0  014791  1015 


0  ++ 


44619.98716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 


D 

SPI 

w 

HT 


0. 

0. 

0 

0 


F107  0. 
FBAR  0. 
AKP  0. 
TVE  0. 
UJD  0. 


ETU 

HO 

GO 

ES 

XI 

RND 


0. 

0. 

0. 

0. 

0. 

0. 


ESD  0. 
AO  0. 
ADOT  0. 


15.71356734  .000783912 
000000000-0  0 
01.06000000  -701.8652  E 

X  xxxxxx.xxxxxxxxx 
0  000001.000000000 
0  000000.000010000 
0  000000.100000000 
0  000200.000000000 
0  000005.000000000 
0  000100.000000000 
0  000050.000000000 
0  000050.000000000 
0  000000.000000000 
0  038422.000000000 
0  043223.000000000 
0  000000.000000000 
0  000000.000000000 
0  000000.000000000 
0  000000.010000000 
0  000022.500000000 
0  000000.000010000 
0  000000.000000000 
0  000001.025940000 
0  000000.000000000 


1.00112  -3.81308  -0.23959E-00 
0  0  0 

00 

xxxxxx.xxxxxxxxx 

000003.500000000 
000000.001000000 
000010 .000000000 
003000.000000000 
000200.000000000 
100000.000000000 
000300.000000000 
000300 . 000000000 
000007 .000000000 
039422.000000000 
044223.000000000 
000179.000000000 
000360 . 000000000 
000360 . 000000000 
000000.100000000 
000157.500000000 
000000 . 000000000 
000000 . 000000000 
000001 . 168470000 
000000.000000000 


++ 

++ 

++ 


1000,  0.1 
2000,  0.3 
3000,  0.5 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


X  X 
X  ONLY  THE  ABOVE  IS  ACTUAL  INPUT  FILE  DATA.  AN  INPUT  ELEMENT  X 
X  DESCRIPTION  IS  PRESENTED  NEXT  TO  AID  THE  USER  IN  PLACING  VALUES  X 
X  x 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT  FILE  FOR  PROPELLANT  LONGEVITY  MODEL  (VARIABLE  ORGANIZATION) 

IETR1 ,  IETR2,  NOA,  PSENS,  IET ,  PSIZE 
CD,  A,  D,  SPI,  W,  HT, 

F107 ,  FBAR,  AKP,  TVE, 

1  SKIP  LINE  ONE  NAVSPASUR  FIVE  CARD  DATA 

2  8X,  UJD,  IX,  ETU,  IX,  HO,  IX,  GO,  IX,  ES,  IX,  XI 

3  20X,  RND,  19X,  ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  3X,  AO,  IX,  ADOT 

IRV(l),  DA 1 1 )  | 

IRVC2),  DAC2)  I  NOTE:  WARNING  OMIT  ALL  ORBIT  ADJUST  DATA  IF 
IRVC3) ,  DA ( 3 )  !  NO  ORBIT  ADJUSTS  ARE  SCHEDULED  CNOA=0) 


IRV(K) ,  DACK)  I 

SENSITIVITY  ANALYSIS  TABLE 

CHAR  SNSAR  SNSMIN  SNSMAX 

A4  F3.1  F16.8  F16.8 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =  1.0000 

F107  (DECIMETRIC  SOLAR  FL^'X)  *  100.00 

FBAR  (90  DAY  AVG.  F107)  =  100.00 

AKP  (GEOMAGNETIC  INDEX)  *  2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =  43222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =  2.20  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 
SPI  (MOTOR  SPECIFIC  IMPULSE)  =  230.00  SECONDS 

H  (INITIAL  FUEL  MASS)  =  10.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  =  9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  44619.98716775  MEAN  JULIAN  DAYS 
ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  *  95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 
XI  (MEAN  INCLINATION)  *  95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =  0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  »  -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =  1.06000  MEAN  EARTH  RADII 

ADOT  (ScMI-MAJOR  AXIS  DECAY  RATE)  »  -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION  5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS  0 


SENSITIVITY  ANALYSIS  FOR 
LOW  EARTH  ORBIT  SATELLITES 
PROPELLENT  LONGEVITY 


TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO 
MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN. 

THE  GRANULARITY  IS  1/  25. 

ALL  OTHER  INPUT  PARAMETERS  ARE  RETURNED  TO  THEIR 
INITIAL  VALUES  AT  EACH  INCREMENT. 


INPUT  PARAMETER  MINIMUM  MAXIMUM  INCREMENT 

CD  1.000000000  3.500000000  0.100000000 


SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS 


REV  NUMBER 

TIMECDAYS) 

FRACTION  OF  SCAN 

495 

31.58 

1.000000000 

448 

28.57 

1.100000000 

409 

26.08 

1.200000000 

376 

23.97 

1.300000000 

348 

22.18 

1.400000000 

324 

20.65 

1.500000000 

303 

19.30 

1.600000000 

285 

18.15 

1 .700000000 

268 

17.07 

1.800000000 

254 

16.17 

1.900000000 

241 

15.34 

2.000000000 

229 

14.57 

2.100000000 

218 

13.87 

2.200000000 

208 

13.23 

2.300000000 

200 

12.72 

2.400000000 

191 

12.14 

2.500000000 

184 

11.70 

2.600000000 

177 

11.25 

2.700000000 

173 

10.30 

2.800000000 

154 

10.42 

2.900000000 

159 

10.10 

3.000000000 

154 

9.78 

3.100000000 

149 

9.46 

3.200000000 

144 

9.14 

3.300000000 

140 

8.88 

j . 400000000 

136 

8.63 

3.500000000 

05000,001,0,1,0,025 

2.2,  .0000500,  1.0,  230. ,001. ,09100. 

100.,  100.,  2.0  ,43222.7382 

1  10947  U  90070  78  060  A  0  1.0  014791  1015  0  0  ++ 

2  10947  44619.98716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 

3  10947  15.71356734  .000733912  1.00112  -3.81308  -0.23959E-00  0  ++ 

4  10947  000000000-0  0000  ++ 

5  10947  01.06000000  -701.8652  E-00  ++ 

cccc  x.x  xxxxxx.xxxxxxxxx  xxxxxx.xxxxxxxxx 

CD  0.0  000001.000000000  000003.500000000 

A  0.0  000000.000010000  000000.001000000 

D  0.0  000000.100000000  000010.000000000 

SPI  0.0  000200.000000000  003000.000000000  ' 

W  0.0  000005.000000000  000200.000000000 

WT  0.0  000100.000000000  100000.000000000 

F107  0.0  000050.000000000  000300.000000000 

FBAR  0.0  000050.000000000  000300.000000000 

AKP  0.0  000000.000000000  000007.000000000 

TVE  0.0  038422.000000000  039422.000000000 

UJD  0.0  043223.000000000  044223.000000000 

ETU  0.0  000000.000000000  000179.000000000 

HO  0.0  000000.000000000  000360.000000000 

GO  0.0  000000.000000000  000360.000000000 

ES  0.0  000000.010000000  000000.100000000 

XI  0.0  000022.500000000  000157.500000000 

RND  0.0  000000.000010000  000000.000000000 

ESD  0.0  000000.000000000  000000.000000000 

AO  1.0  000001.025940000  000001.168470000 

ADOT  0.0  000000. 000000000  000000.000000000 

1000,  0.1 
2000,  0.3 
3000,  0.5 


XXXXXXXXXXXXXXXXXXXXXXXXHXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

X  ONLY  THE  ABOVE  IS  ACTUAL  INPUT  FILE  DATA.  AN  INPUT  ELEMENT 
x  DESCRIPTION  IS  PRESENTED  NEXT  TO  AID  THE  USER  IN  PLACING  VALUES 
x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT  FILE  FOR  PROPELLANT  LONGEVITY  MODEL  (VARIABLE  ORGANIZATION) 

IETR1 ,  IETR2,  NOA,  PSENS,  IET,  PSIZE 
CD,  A,  D,  SPI,  W,  WT, 

F107 ,  FBAR.  AKP,  TVE, 

1  SKIP  LINE  ONE  NAVSPASUR  FIVE  CARD  DATA 

2  8X,  UJD,  IX,  ETU,  IX,  HO,  IX,  GO,  IX,  ES,  IX,  XI 

3  20X,  RND,  19X,  ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  8X,  AO,  IX,  ADOT 
IRV<1),  DA  Cl)  | 

IRVC2),  DAC2)  |  NOTE:  WARNING  OMIT  ALL  ORBIT  ADJUST  DATA  IF 
IRV<3),  DA( 3 )  |  NO  ORBIT  ADJUSTS  ARE  SCHEDULED  CNOA  =  0) 


IRV(K) ,  DA(K)  | 

SENSITIVITY  ANALYSIS  TABLE 

CHAR  SNSAR  SNSMIN  SNSMAX 

A4  F3.1  F16.8  F16.8 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


*  *  *  * 


INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =  1.0000 

F107  (DECIMETRIC  SOLAR  FLUX)  =  100.00 

FEAR  (90  DAY  AVG.  F107)  =  100.00 

AKP  (GEOMAGNETIC  INDEX)  =  2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =  43222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =  2.20  UNITLESS 

A  (CROS SECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 
SPI  (MOTOR  SPECIFIC  IMPULSE)  =  230.00  SECONDS 

H  (INITIAL  FUEL  MASS)  =■  1.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  *  9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  *  44619.98716775  MEAN  JULIAN  DAYS 
ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  *  95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 
XI  (MEAN  INCLINATION)  =  95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =  0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =  -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =  1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =  -701.365200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER. ITERATION  5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS  0 


l 
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TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO 
MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN. 

THE  GRANULARITY  IS  1/  25. 

ALL  OTHER  INPUT  PARAMETERS  ARE  RETURNED  TO  THEIR 
INITIAL  VALUES  AT  EACH  INCREMENT. 


INPUT  PARAMETER  MINIMUM  MAXIMUM 


INCREMENT 


AO 


1.025940000  1.168470000 


0.005701200 


SENSITIVITY  SCAN  RESULTS 


JMBER 

TIME(DAYS) 

1 

0.00 

1 

0.00 

1 

0.00 

3 

0.12 

6 

0.31 

.12 

0.70 

22 

1.34 

40 

2.51 

69 

4.42 

115 

7.47 

186 

12.21 

298 

19.76 

473 

31.65 

747 

50.42 

1149 

78.19 

1652 

113.33 

2033 

140.55 

2270 

158.15 

2514 

176.50 

2942 

208.13 

3349 

238.73 

3643 

262.00 

4026 

291.32 

4373 

318.80 

4555 

334.53 

4712 

348.62 

FOR  ABOVE  SPECIFICS 

FRACTION  OF  SCAN 
1.025940000 
1.031641200 
1.037342400 
1.043043600 
1.048744800 
1.054446000 
1.060147200 
1.065848400 
1.071549600 
1.077250800 
1.082952000 
1.088653200 
1.094354400 
1.100055600 
1.105756300 
1.111453000 
1.117159200 
1.122860400 
1.1285616CQ 
1 .154262300 
1.139964000 
1.145665200 
1 . 151366400 
1.157067600 
1.162768300 
1.163470000 
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05000, 
2.2,  . 
100., 

1  10947 

2  10947 

3  10947 

4  10947 

5  10947 
CCCC  X. 
CD  0 . 


001,0,1,0,025 

0000500,  1.0,  230. ,100. ,09100 
100.,  2.0  ,43222.7382 
U  90070  78  060  A  0 


1.0  014791  1015 


0  ++ 


44619.93716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 


15.71356734 

000000000-0 

01.06000000 


.000783912 
0 


A 
D 

SPI 

w 

WT 


0. 
0. 
0  . 
0  . 
0. 


F107  1. 
FBAR  0. 
AKP  0. 
TVE  0  . 
UJD  0. 


ETU 

HO 

GO 

ES 

XI 

RND 


0  . 
0. 
0. 
0. 
0. 
0. 


ESD  0. 
AO  0. 
ADOT  0. 


X  XXXXXX.XXXXXXXXX 
0  000001.000000000 
o  000000.000010000 
0  000000.100000000 
0  000200.000000000 
0  000005.000000000 
0  000100.000000000 
0  000050.000000000 
0  000050.000000000 
0  000000.000000000 
0  038422. 000000000 
0  043223.000000000 
0  000000.000000000 
0  000000.000000000 
0  000000.000000000 
0  000000.010000000 
0  000022.500000000 
0  000000.000010000 
0  000000.000000000 
0  000001.025940000 
0  000000.000000000 


■701.8652  E-00 


1.00112  -3.81308  -0.23959E-00 
0  0  0 


++ 

++ 

++ 


XXXXXX.XXXXXXXXX 
000003. 50000QQ00 
000000 . 001000QQ0 
000010. 000000000 
003000.000000000 
000200 . 000000000 
,  000000000 
.000000000 
000000000 
.  00000000  0 
039.422 .000000000 
044223.000000000 
.000000000 
.  000000000 
.  000000000 
.100000000 
000157.500000000 
000000.000000000 
000000 . 000000000 
000001.168470000 
000000.000000000 


100000 . 
000300  . 
000300  . 
000007  , 


000179  . 
000360 . 
000360 . 
000000 . 


1000,  0.1 
2000,  0.3 
3000,  0.5 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X  X 
X  ONLY  THE  ABOVE  IS  ACTUAL  INPUT  FILE  DATA.  AN  INPUT  ELEMENT  X 
X  DESCRIPTION  IS  PRESENTED  NEXT  TO  AID  THE  USER  IN  PLACING  VALUES  X 
X  X 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


INPUT  FILE  FOR  PROPELLANT  LONGEVITY  MODEL  (VARIABLE  ORGANIZATION) 

IETR1 ,  IETR2.  NOA,  PSENS,  IET ,  PSIZE 
CD,  A,  D,  SPI,  W,  WT, 

FI 07 ,  FBAR,  AKP,  TVE, 

1  SKIP  LINE  ONE  NAVSPASUR  FIVE  CARD  DATA 

2  3X,  UJD,  IX,  ETU,  IX,  HO,  IX,  GO,  IX,  ES,  IX,  XI 

3  20X,  RND,  19X,  ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  8X,  A0,  IX,  ADOT 

IRVC1),  DA( 1 )  | 

IRVC2),  D A ( 2 5  I  NOTE i  WARNING  OMIT  ALL  ORBIT  ADJUST  DATA  IF 
IP.VC3),  DAC3)  I  NO  ORBIT  ADJUSTS  ARE  SCHEDULED  (N0A  =  0) 


IRVCK) ,  DACK)  I 

SENSITIVITY  ANALYSIS  TABLE 

CHAR  SNSAR  SNSMIN  SNSMAX 

A4  F3.1  F16.8  F16.8 

xxxxxxxxxxxxxxxxxxxxxkxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  *  1.0000 

F107  CDECIMETRIC  SOLAR  FLUX)  =  100.00 

FBAR  (90  DAY  AVG .  F107)  =  100.00 

AKP  (GEOMAGNETIC  INDEX)  =  2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =  43222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =  2.50  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 
SPI  (MOTOR  SPECIFIC  IMPULSE)  =  230.00  SECONDS 

H  (INITIAL  FUEL  MASS)  =  10.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  »  9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  44619.98716775  MEAN  JULIAN  DAYS 
ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =  95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 
XI  (MEAN  INCLINATION)  =  95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =  0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =  -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =  1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =  -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION  5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS  0 
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TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO 
MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN. 

THE  GRANULARITY  IS  1/  25. 

ALL  OTHER  INPUT  PARAMETERS  ARE  RETURNED  TO  THEIR 
INITIAL  VALUES  AT  EACH  INCREMENT. 


INPUT  PARAMETER  MINIMUM  MAXIMUM  INCREMENT 

FI  07  50.000000000  300.000000000  10.000000000 


SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS 


REV  NUMBER 

TIME( DAYS ) 

FRACTION  OF  SCAN 

392 

24.99 

50 . 000000000 

324 

20.65 

60.000000000 

276 

17.58 

70 . 000000000 

241 

15.34 

80 . 000000000 

213 

13.55 

90 . 000000000 

191 

12.14 

100.000000000 

174 

11.06 

110.000000000 

159 

10.10 

120 . 000000000 

146 

9.27 

130 . 000000000 

136 

8.63 

140 . 000000000 

127 

8.05 

150.000000000 

119 

7.54 

160.000000000 

112 

7.10 

170.000000000 

105 

6.65 

180 . 000000000 

100 

6.33 

190.000000000 

95 

6.01 

200.000000000 

90 

5.69 

210.000000000 

86 

5.43 

220 . 000000000 

3? 

5.13 

230  .  OOOOOOOOG 

79 

4.99 

24Q . 000000000 

76 

4.79 

250.000000000 

73 

4.60 

260 . 000000000 

70 

4 .41 

270 . 000000000 

67 

4.22 

280 . 000000000 

65 

4 . 09 

290 .000000000 

63 

3.96 

300 . 000000000 

INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =  1.0000 

F107  (DECIMETRIC  SOLAR  FLUX)  =  100.00 

FBAR  (90  DAY  AVG.  F107)  *  100.00 

AKP  (GEOMAGNETIC  INDEX)  *  2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  *  43222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =  2.20  UNITLESS 

A  (CROSSECTIONAL  AREA)  *  0.000050000  SQUARE  KILOMETERS 
SPI  (MOTOR  SPECIFIC  IMPULSE)  =  230.00  SECONDS 

N  (INITIAL  FUEL  MASS)  *  150.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  »  9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  ’  44619.98716775  MEAN  JULIAN  DAYS 
ETU  (MEAN  ANOMALY)  =»  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  »  95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  *  200.0000  DEGREES 
XI  (MEAN  INCLINATION)  *  95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =  0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  *  -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  *  1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  *  -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION  5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS  0 


PROPELLANT  LIFE  ESTIMATOR 


CDA3  0.11000E-03  SQ  KM 
SPECIFIC  IMPULSE3  230.00  SEC 
INITIAL  WEIGHT  OF  PROPELLANT3  150.00  KG 


F1073  100.00 

FBAR=  100.00 

TIME  OF  VERNAL  EQUINOX3  0 . 4322273S200000E+05 


MAXIMUM  NUMBER  OF  ORBITAL  REVOLUTIONS’ 


DATA  PRINTED  EVERY  100  ORBITAL  REVOLUTIONS 


KOZAI  MEAN  ELEMENT  SET 

SEMIMAJOR  AXIS3  0 . 675620541 07790E+04  KM 
ECCENTRICITY3  0 . 9656 1 1 329046 15E-02 
INCLINATION3  0 . 949999648496 17 E+02  DEG 
ARGUMENT  OF  PERIGEE3  0 . 19333742463643E+03  DEG 
ASCENDING  NODE3  0 . 95000149668470E+02  DEG 
MEAN  ANOMALY3  0 . 15661152641309E+03  DEG 


REV  NUMBER 

TIMECMJD) 

TIME(DAYS) 

PROPELLANT  REMAINING(KG) 

100 

44626.32 

6.33 

145.28 

200 

44632.71 

12.72 

140.69 

300 

44639.10 

19.11 

136.21 

400 

44645.49 

25.50 

131.83 

500 

44651.88 

31.90 

127.54 

600 

44658.27 

38.29 

123.37 

700 

44664.67 

44.68 

119.33 

800 

44671.06 

51.07 

115.45 

900 

44677.45 

57.46 

111.72 

1000 

44683.84 

63.86 

108 . 07 

1100 

44690.23 

70.25 

104.42 

1200 

44696.63 

76.64 

100.72 

1300 

44703.02 

33.03 

96.94 

1400 

44709.41 

89.42 

93.05 

1500 

44715.80 

95.82 

88.96 

1600 

44722.19 

102.21 

84.60 

1700 

44728.59 

108.60 

79.92 

1800 

44734.93 

114.99 

74.98 

1900 

44741.37 

121.33 

69.93 

2000 

44747.76 

127.77 

64.97 

2100 

44754.15 

134.17 

60.27 

2200 

44760.55 

140 . 56 

55.92 

2300 

44766 , 94 

146.95 

51 . 91 

2400 

44773.33 

153.34 

48.13 

2500 

44779.72 

159.73 

44.59 

2600 

44786.11 

166.13 

41  .  07 

2700 

44792.51 

172.52 

37.58 

2800 

44798.90 

178.91 

34.10 

2900 

44805.29 

185.30 

30.53 

3000 

44811.68 

191.69 

26  .78 

3100 

44818.07 

198 . 09 

22.74 

3200 

44824.47 

204.43 

18.37 

3300 

44830.86 

210.87 

13.66 

3400 

44337.25 

217.26 

8.68 

3500 

44843.64 

223.65 

3 . 58 

3600 

44850.03 

230 . 05 

-1  .  49 

3600 

230.05 

0.000000000 
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